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Introduction 


This  report  is  related  to  the  recovery  of  f  -ivity  anomalies  from  satellite 
altimeter  data.  Several  previous  reports  have  been  prepared  on  the  subject 
by  the  author  (Rapp,  1983a,  1979).  The  emphasis  in  our  previous  studies  has 
been  the  determination  of  mean  gravity  anomalies  in  l*xl*  and  5*  extent.  In 
this  study  our  emphasis  will  be  on  point  gravity  anomalies. 

A  second  difference  of  this  study  from  our  previous  ones  relates  to  the 
data  used.  In  our  previous  studies  we  have  used  only  Geos->3  (Rapp,  1979b)  or 
only  Seasat  (Rapp,  1983a)  altimeter  data.  We  now  have  a  combined 
Geos-3/Seasat  data  set  (Liang,  1983)  that  provides  a  denser  data  set  than 
either  alone.  In  addition  new  adjustments  of  the  Seasat  data  have  taken  place. 
These  adjustments  have  substantially  improved  the  consistency  of  the  altimeter 
data  over  what  was  used  in  our  earlier  studies. 

The  results  presented  here  are  not  the  first  attempt  we  have  made  at  the 
recovery  of  point  gravity  anomalies.  Using  Geos-3  data  we  had  previously 
shown  the  good  agreement  between  ship  measured  gravity  data  and  altimetric 
implied  data.  In  this  report  we  will  look  at  many  more  comparisons  and  carry 
out  an  ocean  wide  estimation  of  gravity  anomalies  on  a  uniform  grid. 


Method 


The  method  of  anomaly  recovery  is  based  on  least  squares  collocation. 
The  anomalies  are  computed  from  the  corrected  sea  surface  height  using  a 
reference  gravity  field  and  a  set  of  covariances.  Specifically  we  have: 


=  Cgh  +  D)-‘  (h  -  hR)  +  AgR  (1) 

~  Cgg  -  Cgh  (Qhh  +  D)  ‘  ^g  (2) 


where 

Ag  predicted  free-air  gravity  anomaly  with  respect  to  the  ellipsoid 

gravity  field; 

h  column  vector  of  the  altimeter  implied  geoid  undulations; 

Cgh  row  vector  containing  the  covariance  (referred  to  the  reference 

field)  between  the  anomaly  being  predicted  and  the  given  geoid 
undulations; 

Chh  square,  symmetric  matrix  containing  the  covariances  (referred  to 

the  reference  field)  between  the  given  geoid  undulations; 

D  error  covariance  matrix  of  the  ’observed’  geoid  undulations 

which  was  taken  to  be  a  diagonal  matrix  whose  elements  corres¬ 
ponded  to  the  square  of  the  standard  deviation  of  the  alti¬ 
meter  measurement; 


1 


expected  mean  square  value  (referred  to  the  reference  field)  of 
of  the  anomaly  predicted; 

predicted  anomaly  standard  deviation; 

gravity  anomaly  and  geoid  luidulation  implied  by  the  reference 
gravity  field. 

Equations  similar  to  (1)  and  (2)  can  be  written  for  the  estimation  of  the  geoid 
undulation.  Our  ultimate  aim  will  be  to  apply  (1)  to  obtain  &g  and  the  geoid 
undulations  at  spherical  grid  interval  with  respect  to  a  defined  reference 
system.  We  will  also  discuss  tests  in  which  30’  mean  gravity  anomalies  were 
computed  in  preliminary  studies  that  were  made  early  in  this  research  effort. 


mg 

^gR.  hR 


The  Ellipsoid  Parameters  Used  in  this  Study 

The  adtimeter  data  we  used  refers  to  the  ellipsoid  of  the  Geodetic  Reference 
System  1980.  The  parameters  of  the  ellipsoid  are: 


a  =  6378137  m 
f  =  1/298.257 


In  the  collocation  solutions  we  intend  to  use,  the  data  vector  (i.e.  the  geoid 
undulations)  should  refer  to  the  best  ellipsoid  parameters.  Currently  the  best 
"a"  value  is  thought  to  be  (Rapp,  1983b). 


a  =  6378136  m 

To  refer  the  existing  sea  surface  heights  to  the  new  ellipsoid  we  must  add  1 
meter  to  such  values.  This  will  be  done  for  all  computations  in  this  paper  so 
that  all  geoid  undulations  given  in  this  paper  will  refer  to  the  following 
ellipsoid: 


a  =  6378136  m 
f  =  1/298.257 

The  appropriate  gravity  formula  roust  also  be  considered  as  a  reference  for 
the  gravity  anomalies.  The  equatorial  gravity  for  GRS80  is: 

ra  =  9.780327  ms"* 

=  978032.7  mgals 


The  equatorial  gravity  for  the  ellipsoid  used  in  our  computations  (to  be 
consistent  with  the  6378136  m  equatorial  radius)  is  therefore  (1  m  corresponds 
to  approximately  0.3  mgal  change  in  gravity): 


7a  =  978033.0  mgals 


Before  the  anomalies  derived  from  the  altimeter  data  are  associated  with  a 
specific  gravity  formula  the  effect  of  the  atmosphere  must  be  considered  to 
obtain  anomalies  that  are  consistent  with  terrestrial  anomaly  estimates  (Rapp, 
1979b).  To  do  this  0.87  mgals  must  be  subtracted  from  the  anomaly  derived 
from  equation  (1). 

All  anomaly  comparisons  that  are  made  in  this  report  between  terrestrial 
and  altimeter  derived  anomalies  have  been  made  after  conversion  to  the  same 
gravity  formula  and  the  application  of  the  atmospheric  correction  to  the 
altimeter  anomalies.  The  anomaly  maps  that  will  be  given  do  not  have  the 
atmospheric  correction  applied. 


The  Sea  Surface  Topography  Model 

The  sea  surface  heights  obtained  from  the  altimeter  data  refer  to  the  ocean 
surface  which  is  not  an  equipotential  surface.  To  obtain  the  geoid,  sea  surface 
topography  must  be  removed  from  the  altimeter  defined  sea  surface  height. 
Actual  estimates  of  sea  surface  topography  are  difficult  to  obtain  on  a  global, 
time  critical  basis.  We  have  therefore  used  the  data  of  Levitus  (1982)  to 
develop  a  sea  surface  topography  data  base  for  our  initial  tests. 

We  specifically  started  with  the  mean  annual  dynamic  topography 
computed  by  Levitus  with  respect  to  the  2250  dbar  surface.  The  data 
consisted  of  33856  l*xl*  values  in  the  ocean  regions.  The  mean  value  of  these 
values  was  2.02  dyn-m.  This  value  is  used  to  define  the  geoid  and 
consequently  was  removed  from  all  of  the  original  values.  The  remaining 
values  in  land  areas  were  assigned  a  value  of  zero  so  that  a  complete  data 
array  consisted  of  64800  values.  For  the  evaluation  of  the  sea  surface 
topography  at  an  arbitrary  point  a  bilinear  interpolation  procedure  was  used. 
The  above  procedure  was  developed  and  implemented  in  a  subroutine  by 
T.  Engelis.  A  discussion  of  sea  surface  topography  using  the  Levitus  data  is 
given  by  Engelis  (1983,  1985). 

It  should  be  noted  that  this  sea  surface  topography  represents  a  long 
term  time  average  ('’  20  years)  and  a  spatial  average  on  the  order  of  700  km. 


Covariance  Models 

The  covariances  that  appear  in  equations  (1)  and  (2)  need  to  be  obtained 
from  the  empirical  data,  from  covariance  models,  or  a  combination  of  both. 
Since  (h  -  hR)  is  the  "observed”  quantity  we  could  directly  compute  the 
covariance  function  from  the  individual  samples  in  a  region.  The  covariance 
function  would  reflect  the  geoid  undulation  signal  information,  altimeter  noise, 
and  the  noise  or  error  contribution  from  the  reference  field. 

The  cross  covariances  (Cagh)  obtained  by  empirical  means 

because  we  do  not  have  the  empirical  gravity  data  to  cross  correlate  with  the 
altimeter  residuals.  The  cross  covariances  must  be  obtained  from  some  models 


that  are  somehow  related  to  the  actual  data  area. 


The  first  theoretical  covariance  function  considered  is  that  developed  by 
Tscherning  and  Rapp  (1974)  on  the  basis  of  an  anomaly  degree  variance  model, 
Cn-  Given  such  a  model  the  covariance  functions  with  respect  to  a  degree  N’ 
reference  model  are  as  follows: 


C(4gp,  Agg) 


)  Cn  3"+^  Pn(cos  i'^) 
n-N+l 


(3) 


C(4gp.  Nq) 


y~  7^  Pn(cos  V'pq; 

n=N+l 


'4 


C(Np. 
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z: 

n=N+l 


_ £□ _ 

fn-1)^ 


Pn(cos  ^pq) 
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For  the  Tscherning/Rapp  model  we  have: 


r  -  A(n-l) 

°  “  (n-2)(n+B) 


(6) 


where  A=425.28  mgaP  and  B=24.  For  points  P  and  Q  located  on  a  mean  sphere 
(as  we  will  assume)  we  have  s=0.999617.  Closed  expressions  for  (3),  (4),  and  (5) 
exist  when  N=3.  For  higher  N,  the  covariances  are  found  by  removal  of  the 
contribution  from  degree  3  to  N.  Subroutine  COVA  (in  Tscherning  et  al.  (1974)) 
can  be  used  for  such  a  computation.  (C(Np,  Nq)  is  shown  in  Figure  1  for  N=20 
and  in  Figure  2  for  N=180.  For  N=20,  the  signal  covariance  is  13.32  m^  and  the 
correlation  length  (i.e.,  s  where  C(s)=Co/2)  is  2t30.  For  N=180,  the  signal 
variance  is  0.22  m^  and  the  correlation  length  is  0t29.  The  smaller  correlation 
length  as  N  increases  is  an  advantage  in  using  higher  degree  reference  fields. 

The  errors  in  the  reference  field  can  be  considerea  knowing  the  error 
degree  variances  of  the  model  (Colombo,  1980).  The  two  models  considered  here 
are  the  GEMS  field  (Lerch  et  al.,  1979)  which  is  complete  to  degree  20,  and  the 
180x180  model  of  Rapp  (1981).  The  error  degree  variances  for  the  GEM9  field 
can  be  computed  from  the  published  accuracy  estimates  of  the  coefficients. 
For  the  Rapp  model  the  published  accuracy  estimates  to  degree  36  have  been 
used.  Beyond  this  a  more  optimistic  estimate  assuming  a  global  anomaly 
accuracy  of  *10  mgals  has  been  used. 


The  error  covariances  are  computed  using  equations  (3,  4,  5)  with  <5cn 
replacing  Cn-  The  error  covariances  for  each  field  is  shown  in  Figure  1  (GEM9) 
and  Figure  2  (0SU81)  for  C(Np,  Nq). 


The  total  covariance  function  for  use  in  equation  (1)  is  the  sum  of  the 
signal  and  noise  covariance  functions.  We  have  (Colombo  1980): 

C(l^)  =  Cp(Tfr)  +  4C(f)  (7) 

where  Cji(f)  is  the  covariance  function  implied  by  the  reference  field  and  AC(f) 
is  that  implied  by  the  potential  coefficient  noise  model.  This  sum  for  C(Np,  Nq) 
is  shown  in  Figures  1  and  2.  The  variance  (with  noise)  for  the  GEMS  field  is 
now  16.3  m‘  with  a  correlation  length  of  about  2t6  which  is  slightly  longer  than 
when  the  noise  was  not  added.  The  variance  (with  noise)  for  the  0SU81  field 
is  now  1.35  m’  and  the  correlation  length  is  2t0  which  is  considerably  larger 
than  the  case  without  noise. 

An  alternative  interpretation  to  equation  (7)  has  been  given  by 
Jekeli  (1985,  private  communication)  who  showed  that  AC  could  be  considered 
part  of  the  noise  matrix  D.  To  show  this  Jekeli  represents  the  residual 
altimeter  observation  as 

h  -  hR  =  h  -  hR  +  £h-hR 

where  the  "  indicates  the  observed  quantity  and  e  is  the  error  in  the 
observations.  We  have: 

cov(h-hR,  h-hR)  =  cov(h-hR,  h-hR)  +  cov(eh-hp,  ch-hR) 

If  we  assume  no  correlation  between  ejj  and  ehp  have: 

cov(£h-hR.  ^h-hfj)  =  cov(ch.  «h)  +  cov(ehR.  «hfj) 

We  thus  have: 

cov(h-hR,  h-hR)  =  cov(h-hR,  h-hR)  +  cov(eh,  ch)  +  cov(ehjj>  'hp) 

This  result  is  the  same  as  (7)  when  the  D  given  in  (1)  is  added. 

The  covariance  model  has  two  primary  parameters  (A,  s)  that  can  easily  be 
changed  in  order  to  fit  an  empirical  covariance  function.  The  variance  of  the 
quantity  is  primarily  dependent  on  A  so  that  knowing  one  we  can  estimctd  the 
otner.  The  second  parameter  s  depends  on  the  ratio  of  the  Bjerhammer  sphere 
radius  to  the  mean  radius  of  the  earth: 

-  iri'  (B, 

One  quantity  that  is  important  in  the  modelling  of  an  empirical  covariance 
function  is  the  correlation  length  (Moritz,  1980,  p.  174),  <*.  f  is  the  value 

where  C(c)“C(0)/2.  We  have  computed  this  correlation  length  for  the 
undulation  covariance  function  for  the  Tschernlng/Rapp  degree  variance  model 


using  different  values  of  Rr  (with  R=6371.0  km)  for  a  degree  20  and  a  degree 
180  reference  field.  These  results  are  shown  in  Table  1. 


Table  1 

Geoid  Undulation  Correlation  Lengths  for  the  Tscherning/Rapp 

Degree  Veu’iance  Model 


Rb  -  Re 

Degree  20  Reference 

Degree  180  Reference 

(kn) 

(degree) 

(degree) 

0.10 

2.30 

0.29 

0.50 

2.30 

0.29 

1.00 

2.31 

0.29 

1.50 

2.31 

0.29 

2.00 

2.32 

0.30 

2.50 

2.32 

0.30 

3.00 

2.33 

0.31 

3.50 

2.33 

0.31 

We  see  that  the  correlation  length  is  insensitive  to  Rb  -  Rfi  <for  the 
differences  considered.  As  the  difference  increases,  the  correlation  length 
increases.  We  would  not  be  able  to  match  well,  with  this  model,  undulation 
covariances  that  have  correlation  lengths  shorter  than  given  in  Table  1. 

Other  degree  variance  models  can  be  used  (Moritz (1980),  Rapp(  1979a))  that 
involve  more  parameters  which  enable  shorter  correlation  lengths  to  be 
modeled. 

Another  covariance  model  that  will  be  useful  is  that  given  by  Jordan 
(1972).  He  proposes  a  third  order  Markov  undulation  model  with  consistent 
covariance  and  cross  coveuriance  for  quantities  of  interest  to  us.  In  our 
application  we  consider  these  models  for  the  residual  observed  quantities,  i.e. 
(x  -  Xref)'  Specifically  we  have: 


6 


where: 


91«j2  is  the  undulation  variance; 

is  the  anomaly  variance; 
r  is  the  PQ  separation 

d  is  the  characteristic  distance 


Ij,Kj  are  modified  Bessel  functions  of  order  j 

The  correlation  distance  {()  of  the  stationary  random  process  is  the 
separation  at  which  the  covariance  has  declined  to  1/e  of  the  variance.  That 
is,  C(f)  =  0.367879C(0).  Each  gravimetric  quantity  will  have  its  own  correlation 
distance  that  is  related  to  the  characteristic  distance: 

CN  =  2.905D  =  undulation  correlation  distance  (11) 

Cg  =  1.3610  =  anomaly  correlation  distance  (12) 

This  Jordan  model  has  two  parameters,  and  D,  that  can  be  determined  to 
yield  fits  of  the  model  to  empirical  covariance  functions.  From  such  a  function 
we  have  immediately  vji}.  We  can  numerically  determine  (  from  the  covariance 
function  and  them  compute  D  from  (11)  when  dealing  with  undulations.  Given 
and  0,  the  anomaly  variance^  can  be  computed  from: 


(13) 


where  y  is  an  average  value  of  gravity.  We  note  for  our  application  of  anomaly 
recovery  the  importance  of  C(Np^Agq).  In  the  Jordan  model  this  cross 
covariance  function  is  quite  sensitive  to  iru  and  D.  The  smaller  D  is,  the  larger 
vg,  and  the  larger  the  cross  covariance  is. 


Empirical  Undulation  Covariance  Estimates 

Given  estimates  of  sea  surface  heights  from  satellite  altimetry  we  can 
compute  at  these  points  the  undulation  implied  by  the  GEM9  field  or  by  the 
0SU81  180x180  field.  We  then  have  a  set  of  residual  heights  for  which  we  now 
wish  to  determine  an  empirical  coveu*iance  function. 

Given  a  quantity  s,  at  location  x,y,  the  isotropic  covariance  function  can 
be  written  as: 


Css(r)  =  E((s(X,Y)  -  s) (s(X+x.Y+y)  -  s) 

"  /x»  +  y» 


(14) 
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Here  s  is  the  mean  (in  the  given  data  set)  of  the  quantity  for  which  the  auto 
correlation  function  is  being  computed.  Equation  (14)  is  written  to  show  the 
removal  of  the  mean  value  of  the  observed  quantity  from  the  data  in  the 
computation.  This  effectively  "centers"  (MoritZf  1980,  p.76)  the  data.  We  will 

look  at  covariance  samples  when  s  is  used  and  when  it  is  set  to  zero.  The 

covariance  functions  computed  from  (14)  will  reflect  the  signal  and  reference 
field.  As  such  it  would  be  comparable  to  C(Y’)  given  in  (7)  plus  the  noise  D. 

Given  a  set  of  altimeter  data  in  a  specific  region,  the  covariance  function 
is  computed  by  taking  the  mean  value  of  the  products  as  shown  in  (14).  The 
distance  to  be  associated  with  a  product  is  used  to  store  and  sum  the  product 
in  a  histogram  unit.  For  example,  one  bin  may  contain  all  products  for  which 
0'.9754^<1*.025.  An  average  Y*  is  calculated  for  each  bin  and  is  thus  associated 
with  eacn  mean  product. 

The  tests  that  we  made  took  data  in  a.)  areas  of  3*x3*  and  b.)  areas  of 

2*x2*.  We  carried  them  out  in  several  different  geographic  areas  with  respect 

to  the  GEM9  field  and  the  OSU81  field  as  reference  gravity  models. 

The  first  test  area  is  that  of  the  New  England  seamounts  in  the  western 
Atlantic.  The  empirical  undulation  covariance  with  respect  to  the  GEM9  and 
0SU81  fields  are  shown  in  Figures  3  and  4  where  the  mean  has  ,  and  has  not, 
been  removed  from  the  data  vector.  For  the  GEM9  case  the  mean  was  -2.2m, 
accounting  for  the  large  difference  in  the  residual  undulation  variance.  For 
the  180  case  the  mean  was  0.15  m  so  that  the  two  curves  are  close  together. 
The  data  for  these  two  figures  was  taken  in  a  2*x2*  block. 

Figure  S  show  a  180  covariance  function  in  the  Louisville  Ridge  area  in  the 
southwest  Pacific.  The  mean  of  the  residual  was  0.09m.  The  data  for  this 
graph  was  selected  in  a  3*x3*  area. 


Figure  3 

Empirical  Undulation  Covariance  Function  in  New  England  Sea  Mount  Area 

(GEMS  reference) 
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Of  more  specific  interest  is  the  relation  of  the  empirical  covariance 
functions  to  the  model  covariance  function  described  earlier.  Figures  6  and  7 
show  empirical  and  model  covariance  functions  based  on  data  in  a  3*x3*  area 
around  the  Puerto  Rico  trench  for  the  GEM9  and  OS U 180  reference  fields 
respectively.  Plotted  in  each  figure  is  the  Tscherning/Rapp  covariance 
function  scaled  so  that  the  sum  of  the  residual  and  noise  variance  would  be 
the  same  as  the  actual  variance,  and  the  Jordan  model  with  oii]*  equal  to  the 
empirical  variance,  and  the  D  value  estimated  from  the  correlation  distance. 
For  the  GEM9  reference,  D  was  0.*17  while  it  was  0.*63  for  the  OSU180 
reference  field.  We  note  the  good  fit  of  the  Jordan  function  for  f  up  to  l.*2 
for  the  180  reference  and  higher  for  the  GEM9  reference. 

Figure  8  is  another  example  using  data  in  the  Louisville  Ridge  area  with 
the  180  reference  field.  The  D  value  in  this  case  is  0.*14. 

These  curves  demonstrate  that  the  Jordan  function  can  be  fitted  to  the 
empirical  function  in  a  simple  manner.  The  Tscherning/Rapp  covariances  show 
a  greater  correlation  than  appears  to  exist  for  the  empirical  data  in  the  areas 
tested.  In  some  cases  it  was  not  possible  to  determine  a  D  value  for  the 
Jordan  model  due  to  the  slow  decay  of  the  covariance  function.  In  some  cases 
0  was  found  to  be  small  (i.e.  <  0.*1).  The  small  D  implies  a  large  og  which  can 
become  unreasonable.  In  such  cases  we  will  set  D  to  0.*1. 


Numerical  Values  of  C(Np.Agq) 

We  now  consider  numerical  values  of  C(Np,AgQ)  based  on  equation  (4)  (with 
both  signal  and  reference  field  noise  contributions)  and  equation  (10).  For 
the  use  of  equation  (4),  the  C(Np,Nq)  and  C(Np,AgQ)  functions  have  been  scaled 
to  give  the  same  variance  for  C(Np,Nq)  as  determined  empirically.  For  the  use 
of  (10),  the  has  been  found  through  examinations  of  the  numerical 

covariance  function.  However  to  show  the  dependence  of  C(Np,AgQ)  on  D,  we 
have  used  a  number  of  different  values.  The  value  of  for  two  different 
areas  was  computed  with  data  referenced  to  the  OSU180  field.  The  first  area 
was  near  the  Puerto  Rico  Trench  where  1.93m*  and  the  D  estimate  was 

0.*133.  The  second  area  was  in  the  North  Sea  where  <7^*  was  0.24m*  and  D  was 
0.*199.  For  these  two  areas  we  have  computed  C(Np,AgQ).  Values  with  V  =  0* 
are  shown  in  Table  1. 


COVflRIPNCE  (XIO  '  )  (m') 

?.0  3.0  4.0  5.0  6.0  7.0 


Empirical  and  Model  Geoid  Undulation  Covariances  in  the  Puerto  Rico  Trench 


Area  (GEMS  reference) 
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COVfiRIRNCE  (m^) 


Figure  7 

Empirical  and  Model  Geoid  Undulation  Covariancea  in  the  Puerto  Rico  Trench 

Area  (OSU  180  reference) 
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Figure  8 

Empirical  and  Model  Geoid  Undulation  Covariances  in  the  Louisville  Ridge  Area 

(OSU  180  reference) 
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Table  1 

C(Np,AgQ),  mgal/m 
OSU180  Reference 
Jordan  Model 


D* 

Puerto  Rico 
Trench  Area 

North  Sea 
Area 

0.120 

94mgal/m 

12]Dgal/ni 

0.133 

86 

0.200 

57 

7 

0.400 

28 

3.5 

0.600 

19 

2.4 

0.800 

14 

1.8 

1.00 

11 

1.4 

The  corresponding  value  of  this  covariance  for  the  Tscherning/Rapp  model 
was  26mgal/ni  for  the  Puerto  Rico  Trench  and  Smgal/m  for  the  North  Sea  Area. 

The  large  values  of  C(Np,4gQ)  for  the  Jordan  model  in  the  Puerto  Rican 
Trench  area  implies  that  an  anomaly  will  be  more  affected  by  a  given 
undulation  residual  than  when  the  Tscherning/Rapp  model  is  used.  In  order  to 
have  the  same  C(Np,AgQ,i(cO)  for  the  Jordan  and  Tscherning/Rapp  models  the  D 
value  would  need  to  be  about  0.*4  which  is  twice  that  found  from  the  actual 
data. 


30*x30*  Mean  Anomaly  Predictions 

We  first  examine  the  predictions  of  30'x30'  mean  free-eiir  anomalies  from  our 
combined  Geos-3/Seasat  data  set.  The  specific  data  used  was  the  Geos-3  data 
that  was  adjusted  to  the  Seasat  data  as  described  by  Liang(1983). 

The  primary  purpose  of  the  tests  was  to  see  the  effect  of  using  the  GEM9 
and  OSU180  reference  fields,  and  to  see  the  effect  of  scaling  the  covariance 
functions  to  the  local  area  instead  of  using  a  uniform  global  covariance 
function.  Unless  otherwise  noted,  the  covariance  function  is  that  of 
Tscherning/Rapp. 


Japan  Trench  Predictions  and  Comparisons 

The  area  of  prediction  was  as  follows:  25<‘4*<29,  142**X*<146*.  This 
was  a  4*x4*  area  containing  64  0.5*x0.5*  blocks.  The  terrestrial  data  used  for 

the  comparison  was  computed  from  lO’xlO’  averages  provided  by  Ganeko 
(private  communication,  1982)  and  described  in  Ganeko(1983).  The  root  mean 
square  free-air  anomaly  in  the  area  is  94mgals. 

The  predictions  were  made  using  the  collocation  procedure  represented  by 
equation  (1).  The  number  of  altimeter  observations  used  was  approximately 
420.  A  single  set  of  such  observations  was  used  to  predict  the  mean  anomalies 
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in  a  area  whose  maximum  size  was  4*x4*.  The  results  of  this  comparison  of 
these  values  to  the  terrestrial  data  is  given  in  Table  2. 


Table  2 

Comparison  of  Terrestrial  and 
Altimeter  Derived  30"x30”  Anomalies 
in  the  Japan  Trench  Area 
(mgals) 


Reference  Field 

Corvariance 

Mean  Diff* 

RMS  Diff 

GEM9 

Global 

0.7 

±15.5 

OSU180 

Global 

-1.0 

16.0 

GEM9 

Local 

1.0 

15.2 

«  (alt- ter) 


The  predicted  accuracy  (from  equation  (2))  of  the  predicted  anomalies 
varied  from  9  to  21  mgals  with  the  GEM9,  local  covariance  prediction.  This 
agrees  quite  well  with  the  RMS  defferences  or  comparisons  with  the  t(!rrestrial 
data.  The  RMS  difference  between  the  two  global  prediction*?  was  .^4.6  mgals  so 
that  the  use  of  two  very  different  reference  fields  had  only  a  minor  effect  on 
the  results. 


Puerto  Rico  Trench  Predictions  and  Compeuriaons 

This  area  of  predictions  was  as  follows:  19^^*<21;  293<X*<<297.  This  was  a 
2*x4*  area  containing  32  0.5*  blocks. 

The  terrestrial  estimates  were  made  using  ship  gravity  data  provided  by 
Watts  (private  communication)  1982).  This  data  was  not  adjusted  by  us  for 
crossover  discrepencies.  The  0.5*  predictions  were  made  using  the  collocation 
procedure  with  GEM9  as  a  reference  field.  The  predicted  accuracies  of  the 
terestrial  0.5*  blocks  from  equation  (2)  varied  from  19  to  46  mgals  with  an  RMS 
of  t  30  mgals.  Improved  estimates  could  be  obtained  by  removing  ship  gravity 
bias  and  drift  terms.  The  root  mean  square  anomaly  in  this  area  is  tl56  mgals. 
The  results  of  the  comparisons  of  various  altimeter  derived  solutions  to  the 
terrestrial  data  are  given  in  Table  3. 


Table  3 

Comparison  of  Terrestrial  and  Altimeter 
Derived  30'x30’  Anomalies  in  the 
Puerto  Rico  Trench  Area 
(mgals) 


Reference  Field 

Covariance 

Mean  Diff* 

RMS  diff 

GEM9 

Global 

2.2 

±12.6 

OSU180 

Global 

6.6 

12.6 

GEM9 

Local 

2.1 

12.4 

t  (alt-ter) 


A  second  prediction  area  was  chosen  in  the  Puerto  Rico  Trench  area  for 
improved  terrestrial  predictions  considering  the  correlation  of  the  free-air 
anomaly  and  depth.  Depth  information  was  obtained  at  5’x5'  intervals. 
Liang(1984»  private  communication)  developed  four  different  techniques  to 
incorporate  the  bathymetry  data  in  the  terrestrial  prediction.  No  improvement 
in  the  comparison  of  the  terrestrial  and  the  altimeter  anomalies  was  seen  when 
the  bathymetry  was  used.  In  fact,  several  of  the  procedures  gave 
substantially  poorer  results.  In  this  second  PRT  region  293<x*s296), 

the  RMS  discrepancy  with  terrestrial  data  was  110.4  mgal. 


Blake  Plateau  Predictions  and  Comparisons 

The  last  30*  test  aurea  was  chosen  for  its  mild  anomaly  field  as  opposed  to 
the  much  rougher  fields  in  the  two  prior  trench  areas.  Four  test  regions  were 
used  with  the  terrestrial  estimates  again  provided  by  Watts  data. 

The  first  region  was  28*4’ *32,  284<X*<287.  The  second  region  was  Zi*4’*28 
284<X*^288.  The  third  region  was  2^*4‘*28,  290<x*^293.  The  fourth  region 
was  2\*4'*21,  290<X*^293.  The  RMS  anomaly  in  this  region  was  on  the  order 
of  135  mgals,  with  the  average  predicted  accuracy  of  the  altimeter  derived 
mean  anomalies  l4  mgals.  The  results  of  the  comparisons  are  shown  in 
Table  4. 

The  tests  described  in  the  above  sections  indicate  the  ability  of  the 
satellite  altimeter  data  to  resolve  30*x30’  anomalies.  The  accuracy  of  the 
determination  in  rough  areas  will  be  on  the  order  of  llO  to  15  mgals  while  in 
smoother  areas  it  will  be  on  the  order  of  15  mgal.  Our  predictions  agree  well 
with  the  terrestrial  data  considering  the  accuracy  of  the  data. 


Table  4 

Comparison  of  Terrestrial  and  Altimeter 
Derived  30’x30’  Anomalies  in  the 
Blake  Plateau  Area 
(mgals) 


Region 

Reference  Field 

Covariance 

Mean  Diff 

RMS  Diff 

GEMS 

Global 

5.4 

±9.1 

GEMS 

Local 

5.4 

8.0 

GEMS 

Local 

2.0 

7.1 

GEMS 

Local 

5.5 

t  IV 

GEMS 

Local 

1.5 

5.5 

Point  Anomaly  Predictions-Preliminary  Tests 


We  now  turn  to  the  prediction  of  gravity  anomalies  at  specific  points  for 
comparison  with  ship  data.  Such  tests  have  bee  performed  earlier  (Rapp, 
1979b,  Figure  14)  but  the  tests  to  be  described  here  will  be  more 
comprehensive. 

The  first  test  was  to  take  a  ship  gravity  profile  as  provideed  by  Watts  and 
calculate  the  altimeter  derived  anomaly  using  a  Geos-3/Seasat  data  base.  The 
profile  was  designated  as  V2901  and  was  used  for  earlier  tests  (Rapp,  1979b). 
The  profile  crosses  the  90*  East  Ridge  in  the  Indian  Ocean.  The  ship  anomalies 
are  shown  as  the  continuous  line  in  Figure  9.  The  profile  is  approximately 
1200km  long  with  the  center  (as  plotted)  at  ♦  =  -6.3*,  X  =  89.4*. 

The  altimeter  derived  anomalies  were  computed  by  T.  Engelis.  The  whole 
line  was  divided  into  9  segments  with  altimeter  data  selected  for  each  segment. 
The  data  was  then  used  in  equation  (1)  to  predict  8  points  on  the  profile 
segment.  For  these  predictions  the  GEM9  reference  field  was  used.  A3*  data 
border  was  used  with  a  maximum  of  420  altimeter  points  selected  for  each 
prediction  segment.  This  data  was  selected  to  be  denser  near  (^1*)  the  line 
segment  and  leas  dense  farther  away.  The  mean  difference  (Watts-altimeter) 
was  2.1  mgals,  with  a  root  mean  square  difference  of  t9.1  mgal.  Although  this 
is  a  good  agreement  between  the  two  data  sets,  it  is  clear  from  Figure  9  that 
there  is  more  high  frequency  information  in  the  ship  data.  If  more  profiles 
were  computed,  power  spectral  analysis  of  the  two  data  sources  could  be 
carried  out  as  was  done  by  Eren  (1980)  when  Geos-3  data  alone  was  used  for 
these  types  of  computations. 
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Figure  9 

Ship  and  Altimeter  (GS/SS)  Gravity  Data  Across  the  90*  East  Ridge 
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We  next  turn  to  a  sequence  of  tests  to  better  define  the  prediction 
techniques  for  point  gravity  anomalies  now  to  be  computed  on  a  grid  of  8*x9*. 
In  these  tests  we  were  specifically  interested  in  a.)  the  8*  size  to  be  used 
with  one  data  selection  and  matrix  inversion  and  b.)  the  covariance  model 
(Tscherning/Rapp  or  Jordan)  to  be  used  in  the  prediction.  The  first  test  area 
was  a  l*xl*  area  in  the  Indian  Ocean  where  good  ship  data  was  available.  The 
data  was  used  to  compute  a  lO'xlO*  or  20’x20'  grid  (8  =  -11.8*,  X  z  88.2*)  of 
anomalies  from  the  Geos-3/  Seasat  data. 


Table  5 

Grid  Prediction  Comparison  in  an 
Indian  Ocean  l*xl*  Block 


8* 

Data  Points 

Border 

Ref  Field 

Mean  Diff 
(mgals) 

RMS  Diff 
(mgals) 

1* 

367 

3* 

GEM9 

9.8 

±14.4 

o* 

404 

1* 

OSU180 

9.8 

14.4 

2* 

245 

1* 

OSU180 

10.2 

15.0 

2* 

344 

0.5* 

OSU180 

11.1 

15.8 

2* 

225 

0.5* 

OSU180 

11.1 

±15.4 

We  see  no  substantial  difference  in  the'  results  presented  in  Table  5. 

Other  tests  were  run  with  8  values  of  3*  and  4*  using  a  0*.5  border  and  the 

OSU  180  reference  field.  The  RMS  difference  of  the  common  prediction 
anomalies  in  the  2*,  404  point  case  was  iO.9  mgal  (8  =  3*,  370  points);  ^5.0 
mgals  (9  =  3*,  292  points);  and  *4.5  mgal  (9  =  3*,  370  points).  We  concluded 

that  increasing  the  value  of  8  had  only  a  marginal  effect  on  the  predicted 

anomalies,  beint^  always  less  than  the  accuracy  of  the  predicted  anomalies. 

A  second  lest  was  run  in  a  region  near  the  Puerto  Rico  Trench.  The 
southern  latitude  of  this  area  was  19.5*  and  the  western  longitude  293.5*. 
Anomalies  in  this  area  ranged  from  40  mgals  to  -360  mgals.  Predictions  were 
made  with  two  reference  fields,  various  prediction  block  sizes  (8),  number  of 
data  points  and  data  borders.  Comparisons  were  then  made  to  the  grid 
estimate  based  on  the  Watts  data.  Some  of  these  tests  are  described  in 
Table  6. 
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Table  6 

Grid  Prediction  Comparisons 
in  Puerto  Rico  Test  Region  One 


e* 

Data  Points 

Border 

Ref  Field 

Mean  Diff 

RMS  Diff 

(■gals) 

(mgals) 

1* 

449 

3* 

GEM9 

-0.5 

±24 

1* 

449 

3* 

OSU180 

-0.4 

±20 

1* 

260 

0.5* 

OSU180 

-3.4 

±16 

2* 

402 

0.5* 

OSII180 

-3.2 

±17 

2* 

271 

0.5* 

OSU180 

-4.9 

-16 

3* 

405 

0.5* 

OSU180 

-6.2 

±14 

3* 

291 

0.5* 

OSU180 

-7.7 

-13 

3* 

411 

0.5* 

OSU180 

-9.6 

±17 

From  these  tests  we  decided  that  the  largest  size  b  ock  in  which  a  grid 
prediction  should  be  made  from  one  data  inversion  was  This  was  d 

the  selection  of  the  OSU180  field  as  the  reference  field,  and  a  0.5*  border  of 
data.  The  border  is  taken  from  the  outer  edge  of  the  limits  of  the  grid  being 
predicted. 

We  now  moved  to  teat  the  prediction  accuracy  in  four  other  areas.  The 
primary  intent  of  this  prediction  was  to  consider  which  covariance  function 
(the  Tscherning/Rapp  or  the  Jordan)  gave  better  predictions.  In  this  test  the 
Tscherning/Rapp  covariances  were  scaled  such  that  the  undulation  variance  of 
the  actual  data  was  the  same  as  the  model  undulation  variance.  The  Jordan 
parameters  were  found  from  the  empirical  covariances  as  described  in  a 
previous  section.  The  minimum  D  value  accepted  was  0.1*  (smaller  values 
simply  implied  an  unreasonably  large  anomaly  variance). 

For  all  the  following  tests  the  data  selection  was  made  so  that  when  using 
the  combined  data  set  an  equal  (approximately)  number  of  Geo8-3  and  Seasat 
observations  were  used.  If  this  was  not  done,  one  would  generally  get  more 
Geos-3  observations  because  of  their  greater  density  in  the  areas  tested. 

The  first  area  in  this  group  test  was  again  the  Puerto  Rico  Trench  area. 
Figure  10  shows  the  location  of  the  available  ship  data,  and  Figure  11  shows 
the  anomaly  field  based  on  ship  data.  The  grid  interval  for  Figure  11  and 
most  other  maps  discussed  in  the  test  section  was  0.25*  (decreasing  this 
spacing  to  0.125*  increased  the  computing  time  by  only  5%  with  some  increase 
in  resolution). 

The  anomaly  and  undulation  predictions  were  carried  out  using  the 
procedures  described  above.  The  anomaly  map,  when  using  the. 
Tscherning/Rapp  covariances,  is  shown  in  Figrue  12.  Figure  13  shows  the 
anomaly  map  with  the  Jordan  covariance  function.  The  geoid  undulation  map  is 
shown  in  Figure  14  for  reference  purpose  based  on  the  Tscherning/Rapp 
model. 
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Figure  11 

Gravity  Anomalies  in  the  Puerto  Rico  Trench  Area  Based  on  Ship  Data 

(C.I.  =  10  mgals) 
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Figure  12 

Gravity  Anomalies  in  the  Puerto  Rico  Trench  Area  from  Altimeter  Data  Using 
the  Tscherning/Rapp  Covariances  (C.I.  =  10  mgals) 


Using 


Visual  comparison  of  Figures  12  and  13  shows  them  to  be  quite  close  in 
form,  except  for  more  apparent  high  frequency  information  in  the  northern 
section  when  the  Jordan  model  was  used.  Comparison  with  Figure  11  (the  ship 
data)  shows  better  agreement  with  the  Tscherning/Rapp  map  in  the  sense  that 
the  loop  seen  in  the  northern  part  of  Figure  13  is  not  apparent  in  Figures  11 
and  12. 

Some  statisticeil  comparisons  between  these  data  sets  are  given  in  Table  7. 

Table  7 

Anomaly  Comparison  in  the 
Puerto  Rico  Trench  for  Two  Covariance  models 

(mgals) 


1  Altimeter  Anomalies 

Ship  Data 

T/R  Covariance  Jordan  Covariance 

Mean  Diff  RMS  Diff  1  Mean  Diff  RMS  Diff 

-7.4  ±21 

-8.7  ±27 

The  result  implies  that  the  Tscherning/Rapp  model  will  give  better  results 
in  the  prediction  of  anomalies  although  the  Jordan  model  gave  a  better  fit  to 
the  empirical  covariance  functions.  The  mean  difference  of  about  -7mgals  may 
bo  related  to  a  base  station  question  or  most  probably  to  the  very  large 
('’-360mgal8)  of  some  anomalies  in  this  area. 

The  geoid  undulation  from  the  two  predictions  were  compared  with  a  mean 
difference  of  1  cm  and  a  RMS  difference  of  il6  cm.  A  few  differences  could 
reach  a  magnitude  of  60  cm.  Another  test  was  done  in  the  area  restricting  the 


southern  latitude  limit  of  the  test  region  to  19.75*.  This  reduces  from  the  data 
the  very  large  anomaly  flucuations  apparent  in  Figures  11,12,  or  13.  In  this 
case  the  RMS  difference  (ship  vs.  altimeter)  was  ^14  mgals  using  the  Jordan 
model  and  111  mgals  using  the  Tscherning/Rapp  model. 


The  second  main  test  area  for  the  point  predictions  was  in  the  Blake 
Plateau  area.  Specifically  the  area  in  which  gridded  predictions  took  place 
using  a  uniform  Geos-3/Seasat  data  set,  or  only  Seasat  data,  was  28.25*^^^31  * ; 
284*<X^286.75*.  Predictions  were  made  using  the  Watts  ship  data  and  the 
altimeter  data  with  the  two  covariance  models.  The  results  of  the  comparisons 
are  given  in  Table  8  for  the  two  cases:  the  uniform  Geos-3/Seasat  data  and 
just  the  Seasat  data. 


Table  8 

Anomaly  Comparisons  in  the  Blake  Plateau 
Region  for  Two  Covsu*iance  Models 
(mgals) 


Altimeter  Anomalies 

Data 

T/R  Covariance 

Jordan  Covariance 

Mean  Diff  RMS  Diff 

Mean  Diff  RMS  Diff 

Geos-3/Seasat 

-6.0  114.3 

Seasat  Only 

-6.5  113.5 

The  results  show  that  the  Tscherning/Rapp  function  gives  better  results 
than  the  Jordan  function.  The  use  of  the  Seasat  data  above  also  yields 
slightly  better  results  than  the  combined  data  set.  The  anomalies  implied  by 
the  Geoa-3/Sea8at  data  are  shown  in  Figure  15  for  the  Tscherning/Bapp 
covariance  function  and  Figure  16  for  the  Jordan  function.  The  maps  are  quite 
similar  but  the  one  with  the  Jordan  function  appears  somewhat  rougher.  The 
maximum  discrepancy  of  the  altimeter  and  ship  data  is  on  the  order  of  30  mgals 
and  occurs  in  an  area  (♦=30*,X=284*)  where  the  altimeter  data  is  quite  sparse. 
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Figure  15 

Gravity  Anomalies  in  the  Blake;  Plateau  Region  From  Altimeter  Data  Using  the 
Tscherning/Rapp  Covariance  (C.I.  =  10  mgals) 


The  agreement  of  the  altimeter  derived  anomalies  and  the  ship  data  could 
perhaps  be  improved  if  the  crossover  discrepencies  in  the  ship  gravity  tracks 
were  removed.  In  a  number  of  cases  these  discrepe  icies  could  reach  20 
mgals. 

The  third  area  in  which  point  anomaly  predictions  were  tested  was  the 
North  Sea.  Our  selection  of  this  area  was  based  on  the  altimeter/anomaly 
study  carried  out  by  Brenneke  and  Lelgemann(1983)  and  the  terrestrial  data 
that  was  available  in  the  form  of  6*xl0'  mean  anomalies  as  described  by  Torge, 
Weber,  and  Wenzel(1983).  The  first  step  here  was  to  compute  a  "point”  0.25* 
data  grid  from  the  Torge  et  al  data.  The  predictions  were  done  regarding  the 
6’xlO'  as  point  values  and  using  a  scaled  covariance  function.  The  predicted 
accuracy  of  the  values  is  on  the  order  of  tl  mgals. 

The  specific  region  in  which  the  predictions  took  place  was  45<^*<60; 
0*<X<6*.  The  predictions  were  done  at  one  time  for  four  3*x3*  areas.  Both 
uniform  Geos-3/Seasat  and  Just  Seasat  data  were  used  with  the  two  covariance 
functions.  Statistics  on  the  comparisons  of  the  576  predicted  anomalies  to  the 
terrestrial  data  set  are  given  in  Table  9. 


Table  9 

Anomaly  Comparisons  in  the  North  Sea  Area 
for  Two  Covariance  Models 
(mgals^) 


Altiaeter  Anosialies 

Data 

T/R  Covariances 

Mean  Diff  RMS  Diff 

Jordan  Covariances 
Mean  Diff  RMS  Diff 

Geos-3/Seasat 
Seasat  Only 

0.9  ±11.5 

1.1  ±  9.4 

1.2  ±13.4 

1.6  ±  9.8 

From  this  table  it  appeeu*s  that  the  scaled  Tscherning/Rapp  covariance 
function  yields  better  results  with  a  given  data  set.  We  also  see  that  the 
Seasat  data  alone  gives  better  results  in  this  area  than  using  the  combined 
data  set.  This  implies  that  the  Geos-3  may  not  be  properly  adjusted  or 
sufficiently  accurate  in  this  somewhat  restricted  region. 

To  examine  this  question  further  we  show  in  Figure  17  the  geoid 
undulations  from  the  Seasat  data  alone  and  from  the  combined  Geos-3  data  set 
using  the  Tscherning/Rapp  coveu*iance  model. 

Although  quite  similar  in  appearance,  the  geoid  undulations  from  the 
combined  Geos-3/Seasat  data  appear  more  irregular  than  the  Seasat  alone  plot. 
Maps  made  just  with  the  Geos-3  data  were  more  irregular  and  i-eflect  the 
larger  crossover  discrepency  that  exists  with  the  Geos-3  data.  The  mean 
difference  and  RMS  difference  between  the  undulations  from  Seasat  alone,  and 
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Figure  17 
Geoid  Undulations  in  the 
(Contour  Interval 

North  Sea  Area 
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from  the  combined  data  set  is  -9cm  and  tl9cm.  The  mean  and  RMS 
differencebetween  predictions  using  the  Geos-3/Seasat  data  set  with  the  two 
different  covariance  function  is  1cm  and  t7cm. 


The  gravity  anomalies  for  this  area  are  shown  in  Figure  18  for  the  Seasat 
alone  and  the  combined  Geos-S/Seasat  data  using  the  Tscherning/Rapp 
covariances. 
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Seasat  Only  Combined  Geos-3/Seasat 

Figure  18 

Gravity  Anomalies  in  the  North  Sea  Area 
From  Alimeter  Data 
(Contour  Interval  =  Smgals) 


These  tests  in  the  North  Sea  area  show  that  one  must  treat  the  combined 
altimeter  data  set  with  care  as  the  use  of  the  Seasat  data  alone  gave  better 
results  than  the  combined  set.  Unfortunately  it  seems  to  be  apparent  from  the 
contour  maps  that  undesirable  noise  is  present  in  the  solution.  The  problem 
most  probably  relates  to  the  reduction  and  adjustment  of  the  altimeter  data, 
especially  for  the  short  Geoa-3  arcs  in  a  region  with  only  data  on  one  side 
(that  is,  there  is  no  substantial  altimeter  data  to  the  east  of  the  North  Sea). 

We  now  turn  to  our  last  test  area  which  is  off  the  coast  of  southern 
California.  The  terrestrial  gravity  data  has  been  taken  from  a  5  mgal  contour 
interval  map  prepared  by  Beyer  et.  al.  (1974).  The  specific  limits  of  the  area 
in  which  the  anomalies  were  predicted  was  32 *.00^4^33 *75;  239*^X^241*.  The 

data  from  the  map  was  digitized  on  a  0.125*  grid.  The  terrestrial  map  is  shown 
in  Figure  19.  There  are  rapid  variations  of  the  anomalies  in  several  portions 
of  this  area  with  much  high  frequency  information  present  in  the  data. 

The  altimeter  anomaly  predictions  were  made  on  a  0.125*  grid  to  match  that 
used  in  the  terrestrial  estimation.  The  anomaly  map  based  on  the 
Tscherning/Rapp  covariance  function  with  the  uniform  Geos-3/Seasat  data  is 
shown  in  Figure  20.  It  is  clear  that  much  high  frequency  gravity  information 
is  missing  although  the  general  structure  of  the  field  is  the  same.  The 
comparison  of  the  terrestrial  altimeter  derived  anomalies  is  given  in  Table  10. 


Table  10 

Anomaly  Comparisons  off  the  Coast  of 
California  for  Two  Covariance  Models 
(mgals^) 


Altimeter  Anomalies 


T/R  Covariances 

Jordan  Covariances 

Data 

Mean  Diff  RMS  Diff 

Mean  Diff  RMS  Diff 

Geos-3/Seasat 


-0.6 


il9.7 


-0.2 


±18.8 
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Figure  20 

Gravity  Anomalies  of  the  Coast  of  California  Prom  Altimeter  Data  (C.I 

mgal) 


The  Accuracy  Computations 


The  accuracy  of  the  predictions  is  computed  from  equation  (2).  The  actual 
value  depends  on  the  data  distribution,  data  noise,  and  the  scale  of  magnitude 
of  the  covariances.  We  will  scale  our  covsu'iance  functions  from  area  to  area  to 
match  model  variance  with  actual  variance  (e.g.  for  residual  geoid  undulations). 
These  values  vEU’y  considerably  from  region  to  region  so  that  with  the  same 
data  one  might  have  predicted  accuracies  that  could  range  from  llO  mgals  to 
iso  mgals.  It  is  therefore  difficult  to  interpret  the  accuracy  estimates  from 
equation  (2)  in  a  global  sense  because  of  the  way  in  which  we  scale  our 
covariance  functions.  Within  a  specific  area,  computed  with  the  same 
covariance  function,  the  accuracy  estimates  can  effect  the  relative  accuracy  of 
the  predicted  quantity. 

The  computation  of  the  accuracy  causes  a  substantial  computational  effort 
over  that  if  just  a  prediction  were  done.  To  see  this  let's  write  (1)  and  (2)  in 
the  following  form: 

^8  =  Qgh  •  A  +  Agfj  fl5) 

where ; 

A  =  (Chh  ^  D)"‘  (h  -hfi) 

Given  a  data  set,  the  vector  A  can  be  evaluated  once  before  any  prediction 
actually  takes  place.  The  predictions  then  go  quite  fast  because  each 
prediction  is  made  as  the  product  of  two  vectors. 

Unfortunately,  such  a  modification  of  (2)  is  not  possible  because  the  Cg^ 
vector  depends  on  the  location  of  the  point  to  be  predicted.  In  our  past 
computations  we  have  computed  both  Ag  and  mg  by  writing  (1)  and  (2)  in  the 
form: 


Ag  =  B(h  -  hR)  Agj^ 

(16) 

m^g  =  Cgg  -  BCjjg 

(17) 

®  =  Cgh  +  D)"‘ 

(18) 

B  must  be  computed  for  each  prediction  point  although  the  matrix  inversion 
need  only  be  done  once  for  a  given  data  array. 

If  the  accuracy  estimate  is  not  computed,  the  calculations  of  the  gravity 
anomaly  is  about  eight  times  faster  than  the  calculations  of  both  quantities.  If 
a  vector  processing  computer  were  available  this  accuracy  estimate  penalty 
could  be  susbstantially  reduced. 


It  is  possible  to  obtain  an  approximate  accuracy  estimate  by  considering 
values  rigorously  computed  with  a  specific  covariance  function  and  then 
scaling  this  result  to  another  area  knowing  the  relationship  in  the  variance  of 
the  covariance  functions. 

A  number  of  tests  were  made  computing  the  formal  accuracy  of  the 
predicted  anomalies.  The  average  accuracy  of  a  predicted  anomaly  was  i23mgal 
when  the  residual  anomaly  variance  was  1150.14mgal*.  The  average  accuracy  of 
a  predicted  geoid  undulation  was  10.2m  when  the  residual  undulation  variance 
was  1.60m^.  For  our  prediction  computation  we  computed  the  scaled  accuracy 
estimate  for  each  region  which  was  assigned  to  be  the  same  for  each  point 
predicted  in  the  region.  Specifically  we  used  the  following: 

m(»g)c  =  23[  I 

A  disadvantage  of  this  procedure  is  that  we  are  not  able  to  recognize 
poor  predictions  caused  by  the  lack  of  altimeter  data  in  a  portion  of  a 
region. 

Another  way  will  be  sought  by  examination  of  the  actual  numbers  (or 
count)  of  the  altimeter  in  geographic  cells.  Specific  details  are  found  in 
another  section. 


The  Choice  for  Prediction  Runs 

There  are  a  number  of  decisions  to  be  made  when  choosing  the  final 
prediction  procedures.  Our  tests  described  in  the  previous  section  have  been 
described  to  indicate  which  variables  should  be  selected.  We  choose  the 
following: 

Reference  Field:  OSU  180x180  (1981) 

Covariance  Function:  Scaled  Tscherning/Rapp 
Block  Size;  3*x3*  (for  one  matrix  inversion) 

Data  Numbeer:  Approximately  300  points 
Data  Border:  0.5*  about  block  border 
Grid  Interval:  0.125* 

The  grid  interval  for  the  predictions  made  for  most  of  the  tests  previously 
described  was  0.25*.  However,  decreasing  this  interval  to  0.125*  only  causes 
about  a  5X  increase  in  the  computer  time,  so  that  the  smaller  interval  was 
selected  although  it  may  not  be  justified  on  the  basis  of  the  spectral  content 
of  the  data.  This  will  be  discussed  later  in  the  report. 


The  Satellite  Altimeter  Data  Set 

The  first  global  gravity  field  analysis  that  we  did  with  altimeter  data  was 
with  Geo8-3  data  (Rappi  1979b).  The  second  global  analysis  was  done  with  the 
Seasat  data  (Rapp,  1982,  1983a).  The  gravity  anomaly  recovery  from  the  Seasat 
data  described  in  Rapp  (1983)  was  slightly  flawed  in  parts  of  the  Pacific 
because  some  long  arcs  were  not  properly  adjusted  in  our  crossing  arc 
adjustment.  An  improved  Seasat  adjustment  was  carried  out  in  1982  by  Jaime 
Cruz.  This  second  general  adjustment  eliminated  most  of  the  problem  that 
existed  in  the  first  adjustment.  The  wide  areas  of  large  crossover 
discrepencies  that  existed  in  the  Pacific  now  had  disappeared  (Rapp,  1982, 
Figure  3)  although  there  were  a  few  areas  where  track  dependent  patterns 
were  still  visible.  The  location  of  the  Seasat  data  that  was  available  from  this 
second  crossing  arc  adjustment  is  shown  in  Figure  21. 

In  September  1981  we  received  from  NOAA  a  revised  Geo8-3  data  set  that 
was  acquired  over  a  3.5  year  period.  The  original  data  set  included  10,520 
arcs  starting  from  April  1975  and  extending  to  December  1978.  For  improved 
gravity  field  work  this  data  base  was  adjusted  (using  the  crossing  arc 
technique)  to  the  Seasat  primary  arcs  that  were  initially  used  by  Rowlands 
(1981)  in  our  first  Seasat  adjustment  and  then  by  Cruz  in  the  improved  second 
adjustment.  In  essence  the  Seasat  primary  arcs  were  held  fixed  with  the 
Geos-O  arcs  receiving  bias  and  tilt  correction  based  on  the  crossover 
discrepencies  of  Geos-3  with  the  primary  Seasat  arcs,  a  primary  crossover,  and 
Geos-3  with  itself  (a  local  crossover). 

This  adjustment  was  complicated  by  the  large  amount  of  Geo8-3  data  in 
some  areas  and  by  the  large  number  of  crossovers  in  the  large  time  extent 
data  set.  To  manage  this  problem,  the  world  was  divided  into  six  regions  and 
regional  adjustments  carried  out.  The  regions  were  made  as  large  as  possible 
considering  the  arc  lengths  and  the  data  density.  Although  smaller  regions 
could  give  better  data  fits,  we  are  concerned  that  gravity,  tide,  or  other  types 
of  signal  will  be  absorbed  in  the  parameter  of  the  adjustment.  The 
distribution  of  the  adjusted  arcs  in  the  Geos-3  data  base  are  shown  in  Figure 
22. 

The  complete  details  of  the  combined  Geo8-3/Seasat  adjustment  are  given 
by  Liang(1983).  Before  the  adjustment  of  the  Geo8-3/Sea8at  data  the  RMS 
crossover  discrepancy  was  t2.86m  (Geos-3/Seasat)  and  l3.25m  (Geos-3/Geos-3). 
After  the  adjustments  the  crossover  discrepencies  were  t0.40m  (Geos-3/Seasat) 
and  l0.41  (Geo8-3/Geos-3).  These  values  varied  little  from  region  to  region. 
This  a-posteriori  crossover  discrepancy  of  about  t40  cm  should  be  compared  to 
a  value  of  about  tZB  cm  found  when  just  the  Seasat  data  was  adjusted. 


One  last  modification  was  made  to  our  primary  arc  Seasat  data  base  before 
our  final  dat  set  was  created.  Even  after  our  second  Seasat  adjustment  we 
noticed  53  arcs  where  the  RMS  crossover  discrepency  was  greater  that  i40  but 
less  than  i62cm.  We  felt  that  this  larger  discrepency  was  due  to  a  bias  and 
tilt  adjustment  being  applied  to  arcs  whose  length  was  too  long.  In  order  to 
form  as  consistent  a  data  base  of  altimeter  data  as  possible  we  decided  to 
isolate  those  primary  arcs  on  an  arc  by  eurc  basis.  The  newly  adjusted  data 
then  replaced  the  second  adjustment  sea  surface  heights.  The  typical 
crossover  discrepencies  of  a  "bad"  arc  mth  other  "good"  arcs  was  on  the 
order  of  l36cm  sdter  the  second  adjustment.  At  the  completion  of  the  special 
arc  by  arc  adjustment  the  new  crossover  discrepencies  for  these  arcs  were  on 
the  order  t22cm,  a  considerable  improvement.  A  complete  discussion  of  this 
special  arc  by  arc  adjustment  was  prepared  by  Liang  and  is  in  Appendix  A. 

At  this  point  the  Geos-3  and  Seasat  data  existed  on  separate  tapes  and 
files  representing  the  individual  regional  adjustment  areas.  Although  close  in 
geographic  proximity,  the  regional  adjustment  borders  did  not  match  between 
the  Geos-3  and  Seasat  data.  For  convenience  in  carrying  out  our  prediction 
estimation,  merged  (Geos-3/Seasat)  data  files  were  created  in  convenient 
geographic  blocks.  All  the  altimeter  data  was  referred  to  an  ellipsoid  such 
that  a=6378136m  with  no  sea  surface  topography  corrections  being  made  at  this 
point. 


Sea  Surface  Topography  Considerations  for  Production  Computations 

In  an  earlier  section  we  described  how  a  sea  surface  topography  model  was 
created  for  test  purposes  by  interpolation  from  the  Levitus  estimate  after 
removal  of  a  mean  value.  We  now  consider  the  determination  of  sea  surface 
topography  from  our  Seasat  altimeter  data  and  the  gravity  field  described  by 
the  potential  coefficient  set  of  Rapp  (1981).  A  discussion  of  several  estimations 
of  sea  surface  topography  is  given  in  Appendix  B  prepared  by  Engelis.  In  this 
Appendix  an  SST  data  set  based  on  Levitus’  oceanographic  estimates  is 
described  as  well  as  estimates  based  on  using  Seasat  data  and  the  GEML2 
gravity  field  (Engelis, 1984).  In  addition  a  spherical  harmonic  development  to 
degree  20  of  SST  is  computed  combining  the  Seasat/GEML2  results  to  degree  6 
with  coefficients  from  the  harmonic  analysis  of  the  Levitus  data.  We  have 
compared  the  Levitus  estimates  (from  linear  interpolation  of  values  on  a  l*xl* 
grid)  to  estimates  from  the  spherical  harmonic  expansion  to  degree  20  of  SST, 
and  the  estimates  from  Seasat  minus  the  geoid  undulations  implied  by  the  180 
spherical  harmonic  expansion  of  Rapp  (1981).  The  average  values  were 
computed  in  40*x30*  and  40*x40*  areas.  The  results  eure  shown  in  Figure  23. 

In  some  cases  the  Levitus  values  agree  quite  well  with  those  values 
derived  from  altimetry  and  gravity  field  information.  See,  for  example,  the 
blocks  whose  northwest  corner  are:  0*,160*W;  0*,4*W;  -40*,  169*W;  -40*,120*W; 
-40*,80*W;  -40*,40*W;  -40*,0*;  -40*,40*;  -40*,80*;  -40*, 120*.  There  are  blocks 
where  the  agreement  is  poor:  40*,80*W;  40*,40*W;  70*, 160*’,  40*, 120*; 
40*, 160*.  The  mean  difference  and  standard  deviation  between  A  and  B 
(Levitus)  is  14l42cm.  The  corresponding  vsdues  for  A  and  C  (the  spherical 
harmonic  representation  of  the  SST)  is  17tl5cm.  A  and  C  agrees  better  than  A 
and  B  because  C  was  derived  from  the  Seasat  altimeter  data. 


Figure  23 

Mean  Sea  Surface  Topography  KsLiniales 


It  is  apparent  from  Figure  23  that  the  mean  sea  surface  topography  effects 
can  be  substantial.  The  most  negative  value  is  -131cm  and  the  most  positive 
value  is  120cm  based  on  the  Seasat/Dec  81  computation. 

In  order  to  obtain  the  geoid,  our  computed  sea  surface  heights  should  be 
corrected  by  subtracting  the  sea  surface  topography  defined  by  the  spherical 
harmonic  expansion  given  by  Bngelis  in  Appendix  B.  The  predicted  sea  surface 
represents  a  mean  surface  defined  by  the  Geos-3  and  Seasat  data.  Certain 
short  wavelength  oceanographic  features  (such  as  cold  rings)  will  be  averaged 
out  in  the  process.  However,  globally  the  long  wave  length  signal  will  remain 
in  the  final  estimate  of  the  sea  surface  heights. 

The  gravity  anomalies  derived  from  the  sea  surface  heights  will  be  in  error 
due  to  the  fact  that  the  heights  were  not  corrected  for  sea  surface 
topography.  Since  this  error  should  be  long  wavelength  one  estimate  of  the 
anomaly  error  has  been  obtained  from  the  spherical  harmonic  representation  of 
sea  surface  topography.  A  contour  map  of  the  anomaly  correction  is  shown  in 
Figure  B6  in  the  appendix.  The  maximum  correction  reaches  -1.25  mgal  with 
the  rms  value  of  about  ±0.4  mgal. 

A  series  of  tests  were  run  separately  to  see  the  effect  of  sea  surface 
topography  on  gravity  anomaly  predictions.  These  were  done  using  our 
collocation  procedures,  as  before,  but  now  correcting  the  sea  surface  height 
for  the  sea  surface  topography.  Predictions  were  made  in  three  12*xl2*  areas 
where  l/8*xl/8*  grid  predictions  were  made  using  uncorrected  and  corrected 
sea  surface  heights.  The  three  areas  were  selected  on  the  basis  of  the 
behavior  of  sea  surface  topography  as  seen  in  Figure  B6.  Area  1  (in  the 
eastern  South  Pacific  Ocean)  was  selected  because  of  its  fairly  large  gradient 
of  SST.  Area  2  (in  the  northwest  Paicfic  Oceem)  was  selected  because  of  the 
large  magnitude  (about  130  cm)  of  SST.  Area  3  (in  the  North  Atlantic  Ocean) 
was  selected  because  of  the  nominal  magnitude  (about  30  cm)  of  SST.  The 
results  of  the  comparisons  of  the  two  predicted  anomaly  sets  in  given  in  Table 
11  for  the  9216  points  in  each  area. 

Table  11 

Predicted  Anomaly  Difference'*’  Due  to  Sea 
Surface  Topography  Effects 
(mgals) 


Area 

Geographic  Limits 

Ave  Diff 

RMS  Diff 

Std. 

Max  Diff 

1 

-48*<^<f-60- 

240*^X<252* 

-1.8 

±1.9 

±0.7 

-4.3 

2 

36*<4-s48* 

160*<X<180’ 

2.1 

±2.2 

±0.5 

3.7 

3 

24*<4^36* 

303*<X<315* 

0.5 

±0.6 

±0.4 

1.4 

Results  from  uncorrected  sea  surface  heights  minus  results  from  corrected 
sea  surface  heights. 


From  Table  11  we  see  that  the  rms  (root  mean  square)  difference  is  on  the 
order  of  t2  mgal  which  is  much  below  the  accuracy  of  the  anomaly  prediction. 
However  the  effect  is  systematic  and  can  cause  long  wavelength  error  in  the 
recovered  anomalies.  The  average  difference  results  given  in  Table  11  do  not 
agree  with  the  smaller  values  of  the  corrections  implied  by  Figure  B6  in  the 
Appendix.  The  reason  for  the  discrepancy  is  not  clear.  It  is  clear  that  the 
collocation  solution  can  be  very  sensitive  to  sea  surface  topography  estimates 
depending  on  the  degree  of  the  reference  field.  The  higher  degree  the  field, 
the  better  the  fit  to  the  observation.  The  resulting  residuals  will  be  small  and 
may  be  primarily  due  to  SST  if  not  otherwise  corrected. 

Our  computations  will  be  for  sea  surface  heights  and  gravity  anomalies 
from  uncorrected  sea  surfce  heights  from  altimetry.  To  obtain  the  geoid  from 
our  sea  surface  heights,  we  must  apply  sea  surface  topography  corrections. 
There  is  no  easy  way  to  obtain  corrected  anomalies  using  our  procedures. 
Acceptance  of  the  predicted  values  relies  on  the  t2  mgal  error  estimate 
obtained  from  the  dual  prediction  tests. 


Operational  Data  Selection  and  Production  Runs 


The  testing  of  anomaly  production  procedures  described  earlier  indicated 
that  approximately  300  data  points  were  sufficient  for  the  prediction  in  a  3*x3* 
area.  We  started  our  production  procedures  with  this  approach  with  the  data 
being  selected  from  our  complete  adjusted  altimeter  data  tapes.  There  were 
2,272,307  Seasat  and  3,648,133  Geos-3  values  on  this  tape.  The  data  was  first 
selected  in  the  area  of  interest.  This  set  was  then  reduced  to  approximately 
300  data  points  through  a  selection  of  Geo8-3/Seasat  data  in  the  same 
proportion  that  existed  in  the  complete  data  set  for  the  area.  From  this  data 
we  computed  the  mean  residual  of  (N(altimeter)-N(180))  and  its  standard 
deviation.  We  then  deleted  any  point  where  the  residual  minus  the  mean 
exceeded  2.5  times  the  standard  deviation.  This  process  would  delete  a  few 
bad  points  that  still  remained  in  the  data  set.  Most  of  these  points  were  from 
the  Seasat  data. 


The  results  of  the  above  process  were  not  satisfactory  because  some  track 
error  that  was  still  present  in  our  data  gave  very  clear  signatures  in  the 
anomaly  field  with  much  less  signature  in  the  sea  surface  height  maps.  An 
example  of  this  is  seen  in  Figure  24.  Systematic  anomaly  trends  along  the 
altimeter  arc  are  clearly  apparent.  The  track  causing  this  error  was  identified. 
The  predictions  were  then  repeated  and  the  revised  values  used  in  our  final 
set. 


This  pattern  structure  occured  in  a  number  of  areas,  but  mainly  in  the 
Pacific  Ocean.  In  order  to  reduce  this  problem  a  different  data  selection 
method  was  used.  Specifically  a  "thinned"  data  set  was  created  from  the 
complete  altimeter  file.  The  basis  of  the  thinned  data  base  was  the  selection  of 
one  altimeter  measurement  in  every  lO'xlO’  area  in  the  oceans  where  the 
altimeter  data  was  available.  When  there  was  more  than  one  data  point  in  the 
lO’xlO’  area,  the  sea  surface  heights  were  meaned  and  the  single  point  closest 
to  the  mean  was  selected  for  retention  in  the  thinned  data  file.  The  retention 
was  independent  of  the  altimeter  type  (i.e.  Seasat  or  Geos-3).  The  purpose  of 
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Figure  24 

Anomaly  Contour  Map  Showing  Structures  Caused  by  Unremoved  Track  Error 

(Contour  Interval  is  10  mgal) 
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the  thinning  is  to  remove  (as  much  as  possible)  data  that  would  be  associated 
in  a  sequential  way,  with  a  specific  track. 

Our  final  production  predictions  used  the  thinned  data  tape  that  contained 
522,259  Seasat  and  553,044  Geo8-3  measurements.  The  maximum  number  of 
points  available  for  the  predictions  in  a  3*x3*  block  would  be  576.  Where  the 
number  exceeded  300,  the  actual  points  used  were  reduced  to  300.  The 
standard  deviation  check  was  then  applied  to  the  selected  data  with  "bad” 
points  being  removed.  Where  the  altimeter  data  was  sparse  the  number  of 
available  points  could  be  much  less  than  300.  In  no  case  were  predictions 
made  if  there  were  fewer  than  75  data  points  within  the  3*x3*  prediction  block 
and  it’s  0*.5  borders. 

The  actual  predictions  were  carried  out  in  a  systematic  way.  Generally  a 
single  run  predicted  the  anomalies  and  sea  surface  heights  in  16  3*x3*  blocks. 
The  execution  time  for  a  single  run  depends  on  the  data  available,  but  a 
typical  time  was  3.5  minutes  on  an  Amdahl  470  V/8.  Predictions  near  the  coast 
lines  required  different  size  prediction  blocks.  In  coastal  areas  predictions 
would  also  be  made  in  land  areas  near  the  coast.  Such  values  are  not  to  be 
considered  reliable.  In  those  areas  in  which  sufficient  altimeter  data  was  not 
available  the  anomaly  and  sea  surface  height  were  taken  from  the  Dec  1981 
potential  coefficient  set  to  degree  180. 

At  the  completion  of  each  run  an  anomaly  and  sea  surface  height  contour 
map  were  automatically  prepared.  These  maps  were  checked  for  track  error 
signal  and  bad  data  point  effects.  The  maps  were  also  visually  compared  to 
the  GEBCO  bathymetric  charts,  and  to  anomaly  maps  available  to  us  in  the 
ocean  regions.  These  comparisons  were  made  to  verify  unusual  features  or  to 
locate  possible  substantial  errors  in  our  predictions.  In  a  number  of  cases  we 
found  signals  in  the  predictions  that  were  caused  by  track  er  'or  or  bad  data 
points.  Most  of  the  errors  were  found  from  the  examinations  of  the  anomaly 
maps  and  not  the  sea  surface  height  maps.  This  is  because  t  \e  effect  of  sea 
surface  height  errors  is  greatly  magnified  in  going  to  the  anon  aly  domain. 

Our  response  to  the  obvious  errors  was  the  identification  and  removal  of 
bad  data  arcs  or  data  points  with  the  prediction  runs  repeated.  In  some  cases 
we  repeated  our  computation  with  just  "unthinned"  Seasat  data  and  Geos-3 
data.  This  was  done  to  check  the  reliabililty  of  the  determination  of  an 
unusual  feature.  If  the  feature  was  implied  by  both  data  sets,  it  was 
considered  verified.  In  some  cases  it  was  not  possible  to  identify  specific 
tracks  for  deletion.  However,  almost  all  remaining  track  error  was  in  the 
Seasat  data  so  that  the  predictions  in  some  areas  with  clear  track  error  was 
made  only  with  Geos-3  data. 

Some  prediction  problems  were  noticed  in  coastal  areas  and  in  areas  below 
-60*  to  -65*  in  latitude.  In  such  areas  the  sea  surface  height  maps  were  not 
unreasonable  but  the  anomaly  values  could  change  quite  rapidly.  The  problem 
is  the  southern  most  latitudes  was  primarily  related  to  the  measurements  on 
the  ice  by  Seasat.  Such  values  were  clearly  not  edited  out  of  our  data  set. 
We  do  not  consider  reliable  any  anomalies  below  *=-65*.  In  coastal  regions  the 
anomalies  could  be  effected  by  the  lack  of  data  on  the  coast  side  of  the 
prediction  area.  Predictions  could  be  especially  erratic  in  bays  and  narrow 


waterways.  We  do  not  consider  reliable  any  anomalies  predicted  within  about  50 
km  of  a  coast. 

Our  last  test  for  substantial  track  error  was  made  by  preparing  color 
contour  anomaly  maps  on  a  Tetronix  4115B  terminal.  It  is  much  easier  to 
identify  track  problems  on  a  color  display  than  thru  examinations  of  contour 
maps.  A  few  track  errors  still  present  were  removed  by  the  visual  means. 
Howeverf  track  error  is  still  represented  in  our  final  data  set.  This  error 
shows  up  as  narrow  structures  in  the  anomaly  map.  These  features  are  due  to 
small  (10  mgal)  effects  of  unremoved  track  error.  They  are  most  clearly  seen 
in  the  south  Pacific  where  the  anomaly  field  is  quite  mild  and  the  track  error 
effects  stand  out.  An  example  of  this  is  shown  (unfortunately  not  in  color)  in 
Figure  25. 
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Figure  25 

Image  of  Gravity  Anomalies  Recovered  in  the  South  Pacific 
Showing  Linear  Features  Caused  be  Residual  Track  Error 
Color  Spacing  Interval  was  10  mgals 


The  track  error  remaining  in  the  anomaly  field  can  be  filtered  out  in 
several  different  ways.  One  procedure  tried  with  this  data  was  to  form  a 
weighted  average  of  points  selected  to  be  perpendicular  to  the  direction  of  the 
satellite  track.  A  five  point  average  was  tested  by  Engelis  in  the  same  area  as 
Figure  25  with  the  results  shown  in  Figure  26.  Although  much  of  the  appcurent 
track  error  is  removed  some  remains.  However,  it  is  also  cleeu*  that  high 
frequency  information  is  very  lost.  Although  we  could  have  created  a  point 
data  set  that  had  the  track  error  removed  by  smoothing  we  chose  not  to 
because  of  the  possible  loss  of  high  frequency  data.  Additional  study  needs  to 
be  carried  out  to  develop  procedures  to  filter  out  features  of  a  linear  type 
that  su*e  specifically  in  the  direction  of  the  altimeter  ground  track.  A 
procedure  that  was  applied  to  terrestrial  data  by  Forton  (1984)  may  also  be 
applicable  to  this  problem. 
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Figure  26 

Image  of  Gravity  Anomalies  with  Smoothing 
To  Reduce  Features  Caused  by  Track  Error 


The  Standard  Deviationa  of  the  Predicted  Quantities 


As  discussed  in  an  earlier  section,  we  did  not  compute  specific  accuracy 
estimates  for  each  predicted  point.  Instead  a  uniform  standard  deviation  was 
assigned  using  a  covariance  scaling  technique.  However,  this  value  does  not 
take  into  account  the  distribution  of  the  altimeter  data  used  in  the  prediction. 
To  consider  this  we  used  an  empirical  technique  that  was  based  on  the  number 
(N)  of  altimeter  measurements  in  the  0*.5x0*.5  block  in  which  the  point  fell. 
Specifically  we  used: 


SD(New)  =  SD(Old) 


(21) 


If  N  were  zero,  the  factor  multiplying  SD(old)  was  set  to  be  2.5.  The  location 
of  the  0*5x.0*.5  blocks  in  which  there  is  no  data  is  shown  in  Figure  27.  If  no 
predictions  were  made  in  a  3*x3*  region  because  of  too  few  points  the 
anomalies  and  geoid  undulations  from  the  Dec  81  180  potential  field  were  used. 
In  this  case  the  assigned  standard  deviations  of  99.0  mgals  were  used, 
primarily  to  identify  the  values  as  being  from  a  spherical  harmonic  field.  In 
two  regions  close  to  land  areas  the  anomalies  were  very  erratic  and  unrealistic. 
The  anomalies  in  these  areas  had  their  standard  deviations  set  to  100.0  mgals. 
The  two  regions  were  60*^  ♦  *  54*,  8*<<  X  <12*  ,  and  66*'<  ♦  <57*,  291*<  X 
<297*. 


The  location  of  0*.5x0*.S  blocks  where  at  least  one  point  value  has  a 
standard  deviation  greater  than  30  mgals  is  shown  in  Figure  28.  Many  of  the 
areas  represent  land  areas  in  which  the  predictions  were  automatically  done. 


Statistical  Analysis  of  the  Point  Predictions 

The  total  number  of  predicted  anomalies  and  sea  surface  heights  was 
2,399,833.  However  223,065  of  these  values  had  standard  deviations  greater 
than  30  mgals  which  would  indicate  unreliable  estimates.  (Most  of  these  values 
were  on  land  as  seen  from  Figure  28).  We  therefore  computed  statistics  only 
for  those  values  where  the  standard  deviation  was  <  30  mgals.  In  this  case 
the  mean  value  (uncorrected  for  the  effect  of  the  atmosphere)  was  -0.77  mgals; 
the  root  mean  square  value  was  *  28.6  mgals;  the  maximum  value  was  396  mgals 
and  the  minumum  value  was  -360  mgals.  The  frequency  distribution  of  these 
anomalies  is  shown  in  Table  12. 


Figure  28 

Location  of  OiSxO'S  Blocks  Where  One  or  More  Point  Anomaly  Values  Have  Standard  Deviations  Greater 

Than  30  mgals 


Number 


-400 

to  -350 

5 

-350 

to  -300 

93 

-300 

to  -250 

421 

-250 

to  -200 

1493 

-200 

to  -150 

4044 

-150 

to  -100 

8676 

-100 

to 

-50 

54603 

-50 

to 

0 

1074079 

0 

to 

50 

962535 

50 

to 

100 

60100 

100 

to 

150 

8156 

150 

to 

200 

2069 

200 

to 

250 

408 

250 

to 

300 

61 

300 

to 

350 

17 

350 

to 

400 

8 

The  mean  value  of  the  predicted  standard  deviation  was  *12.2  mgals;  the 
root  mean  square  standard  deviation  was  *12.7  mgals;  the  minimum  value  was 
*10  mgals.  The  maximum  standard  deviation  of  the  original  data  set  was  99 
mgals. 

The  location  of  the  0*.5x0*.5  blocks  in  which  anomalies  exceeded  lOOmgals 
in  absolute  value  (and  standard  deviations  were  less  than  30  mgals)  are  shown 
in  Figure  29.  Clearly  these  vedues  generally  follow  well  known  trench 

structures  in  the  oceans.  A  number  of  these  values  occur  below  -60*  latitude. 
Such  values  are  primarily  due  to  altimeter  measurements  on  ice.  We  were 
unable  to  completely  edit  altimeter  measurements  that  were  contaminated  by  ice 
effects.  The  sea  surface  heights  in  such  areas  did  not  show  clear  problems, 
but  the  anomalies  computed  from  such  data  were  quite  erratic.  Generally  we 
should  view  carefully,  any  anomalies  below  -60*  latitude,  and  completely 
discard  any  anomaly  values  below  -65*  latitude. 

The  maximum  anomaly  that  was  predicated  (i.e.  396  mgals)  occurred  at 
♦=19*.5,  X=204*.5  which  is  just  at  the  eastern  ecfge  of  Hawaii.  Another  group  of 
large  anomalies  (374  mgals)  occured  in  the  western  Pacific  over  the  Ogasawara 
Ridge,  4=28 *.00,  X=142*.00.  The  Ridge  is  located  just  to  the  west  of  the  Izu 
Trench.  Over  the  trench  the  anomalies  are  very  negative.  For  example,  at 
4=28*.00,  Xrl43*.375,  the  anomaly  is  -141  mgals.  This  yields  an  average 
horizontal  gradient  of  3.8  mgals/km. 


Figure  29 

Location  of  0t5x0*5  Blocks  in  Which  a  Point  Anomaly  Exceeded  100  mgals  in  Absolute  Value 


The  most  negative  anomaly  (-361  mgals)  in  the  eighth  degree  data  occurs 
at  ♦=19*.5,  X=293*.50  over  the  Puerto  Rican  Trench.  Considering  the  positive 
anomalies  100  km  to  the  north,  the  horizontal  anomaly  gradient  in  this  case  is 
approximately  3.6  mgals/km. 

For  the  points  in  which  the  anomaly  standard  deviation  was  *30  mgeds  the 
mean  sea  surface  height  was  1.69  m;  the  root  mean  square  value  was  ‘31.1  m. 
The  maximum  sea  surface  height  waa  83.8  m  while  the  minimum  was  -107.1  m. 
The  mean  standard  deviation  was  0.13  m.  The  mean  value  of  1.69  m  implies  that 
the  mean  geoid  undulation  in  the  ocean  areas  must  be  on  the  order  of  -1.65  m 
since  the  mean  sea  surface  topography  i.e.  (f-N)  waa  shown  earlier  (see  Figure 
23)  to  be  approximately  4  cm. 


Selected  Anomaly  and  Sea  Surface  Height  Maps 


The  predicted  anomalies  and  sea  surface  heights  provide  a  uniform  data  set 
for  almost  all  ocean  areas.  The  anomalies  have  a  predicted  accuracy  on  the 
order  of  ‘13  mgals.  The  accuracy  deteriorates  near  coastal  regions  and  in  the 
south  polar  regions.  In  most  other  areas  the  data  provides  a  large  amount  of 
detailed  gravity  field  data.  To  demonstrate  this,  we  selected  a  number  of  areas 
in  which  to  make  contour  maps  from  the  0’.125  data. 


Figure  30  shows  a  perspective  view  of  the  relative  sea  surface  in  the  area 
of  the  Magellan  Seamounts  near  the  Bonin  and  Mariana  Trenches.  The  value 
shown  is  the  predicted  sea  surface  height  minus  the  undulation  implied  by  the 
GEMS  potential  coefficients  to  degree  20.  The  removal  of  the  moat  systematic 
part  of  the  sea  surface  clearly  reveals  the  signatures  of  the  sea  mounts  and 
trenches  in  this  area.  Figure  31  shows  the  gravity  anomalies  in  the  same  area. 
The  gravity  anomalies  shown  in  Figure  31  range  from  -256  mgals  to  185  mgals. 
The  most  positive  anomaly  occurs  at  ♦=15*.5,  X=155*.2  which  is  at  the  center  of 
a  large  seamount  in  the  East  Mariana  Basin.  It  is  also  clear  that  the  anomaly 
gradients  can  be  quite  large  reaching  80  mgals  in  13  km. 


A  second  area  that  is  of  interest  is  near  the  Aleutian  Trench.  The  sea 
surface  heights,  in  perspective  view  are  shown  in  Figure  32.  This  view  clearly 
shows  the  trench,  Bowers  Ridge  and  Basin,  the  Shirshov  Ridge  and  the  start  of 
the  Emperor  Seamount  chain.  The  corresponding  gravity  anomaly  map  is  shown 
in  Figure  33. 


Another  area  showing  large  anomaly  variations  is  the  area  of  the  Sula 
Islands  in  the  Molocca  Sea.  The  anomalies  are  shown  in  Figure  34.  The 
predicted  anomalies  range  from  -256  mgals  to  229  mgals  in  a  distance  of 
approximately  170  km  or  an  anomaly  gradient  of  2.8  mgal/km.  These  large 
variations  may  be  compared  to  similar  large  changes  on  the  terrestrial  map  by 
Watts  and  Bodine  (1978). 
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Figure  31 

Gravity  Anomalies  in  the  Area  of  the  Mariana  Trench  and  Magellan  Seamounts 

(Contour  Interval  is  10  mgala) 
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Another  area  showing  significant  seamount  signature  is  in  the  south  west 
Indian  Ocean.  A  perspective  map  showing  several  sea  mounts  including  the 
Lena  and  Ob  seamount  is  shown  in  Figure  35  where  the  sea  surface  height 
minus  the  geoid  undulation  from  the  GEMS  potential  coefficients  is  displayed. 
The  view  looks  in  a  south  east  direction. 

Another  area  that  shows  the  detail  in  the  predicted  data  is  that  related  to 
the  Mendocino  Fracture  Zone  in  the  north  east  Pacific  Ocean.  The  signature  of 
this  structure  is  clear  in  Figure  36  which  is  a  perspective  view  of  the  sea 
surface  heights  minus  the  undulations  computed  from  the  GEMS  potential 
coefficients.  The  step  structure  across  the  Mendocino  Fracture  Zone  is  clear 
as  well  as  a  smedler  step  across  the  Murray  Fracture  2^ne. 


Lena  Seamount  Signature 


Figure  35 

Perspective  View  of  N-N(GEM9)  in  the  Area  of  the  Lena  and  Ob  Seamounts 

(Units  are  meters) 


Figure  36 

Perspective  View  of  N-N(GEM9)  in  the  Area  of  the  Mendocino  Fracture  Zone 


LATITUDE 


For  our  last  group  of  graphical  displays  of  the  0*.125  data  set  we  look  at 
anomaly  signatures  in  three  relatively  small  areas.  The  first  area  is  in  the 
north  Pacific  in  the  area  of  Johnston  Island  which  is  located  at  ♦=16*.6,  \=190t1. 
The  gravity  anomalies  and  sea  surface  heights  in  the  vicinity  of  this  island  are 
shown  in  Figure  37. 
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Figure  37 

Gravity  Anomalies  and  Sea  Surface  Heights 
In  the  Vicinity  of  Johnston  Island 
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The  last  area  to  be  considered  is  in  the  North  Atlantic  over  a  part  of  the 
New  England  Seamounts.  With  our  original  predictions  the  gravity  anomalies 
over  a  portion  of  the  sea  mounts  are  shown  in  Figure  39. 
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Figure  39 

Gravity  Anomalies  in  New  England  Seamount  Area 
Based  on  Production  Data  Selection  Method 


From  bathymetry  and  previous  studies  (Rapp,  1984)  we  know  that  there  is  a 
seamount  at  #=38*. 8  and  X=296*.2.  It  is  clear  from  Figure  38  that  no  significant 
seamount  signature  is  seen.  Examinations  of  the  altimeter  data  used  in  the 
prediction  run  showed  no  points  selected  in  the  immediate  vicinity  of  the 
known  seamount.  A  second  prediction  run  was  made  using  500  data  points  in  a 
2*x2*  area  surrounding  the  seamount.  The  anomaly  map  is  shown  in  Figure 
40. 
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Figure  40 

Gravity  Anomalies  in  the  New  England  Seamount  Area 
Using  500  Altimeter  Points  in  a  2*x2*  Area  Plus  Border 


We  now  start  to  see  the  seamout  signature  in  the  expected  location  with  a 
maximum  gravity  anomaly  of  about  23  mgals.  This  magnitude  was  less  than 
expected  so  that  a  third  prediction  was  made  using  300  data  points  in  and 
around  the  l*xl*  centered  about  the  seamount.  The  results  of  this  prediction 
are  shown  in  Figure  41.  Now  a  very  strong  signature  is  seen  that  reaches  62 
mgals.  In  fact  the  180  reference  field  is  negative  (about-20  mgals)  in  the  area 
so  that  the  correction  to  the  reference  field  reaches  85  mgals. 

The  point  of  the  above  discussion  is  to  show  how  we  can  miss  the  signal 
association  with  a  high  frequency  phenomenon  using  the  data  selection  and 
prediction  procedures  that  we  have  used.  Based  on  a  general  analysis  we  feel 
that  much  detail  has  been  detected.  However  based  on  this  example  we  know 
some  could  have  been  missed.  The  data  stored  in  our  final  data  tape  was  that 
associated  with  Figure  39.  This  was  done  for  consistency  of  the  predictions 
with  the  surrounding  regions. 

The  primary  intent  of  this  graphical  presentation  of  the  results  has  been 
to  indicate  the  detail  that  is  present  in  the  0.125  data  set.  This  detail  is  seen 
more  readily  in  the  gravity  anomaly  field  than  in  the  sea  surface  height  field. 
Howeveri  at  this  point  we  have  no  quantitative  statement  of  the  information 
content  or  spectral  limit  of  the  data.  We  turn  to  this  in  the  next  section. 


Spectral  Content  of  the  0*. 125x0*. 125  Predicted  Valuea 


The  0*.125  grid  interval  was  originally  chosen  so  that  no  information  would 
be  lost  in  the  prediction/gridding  process.  In  principle  we  wanted  to  retain  as 
much  high  frequency  information  in  our  final  data  set  that  is  consistent  with 
the  original  data  and  our  computational  resources.  We  now  must  see  if  our 
gridding  is  in  fact  a  reasonable  one,  or  should  we  have  chosen  a  larger 
interval  such  as  0*.25. 

The  original  data  was  given  at  one  second  intervals  implying  along  track 
spacing  of  7  km.  Brammer  and  Sailor  (1980)  estimated  the  resolution  limit  of 
the  SEASAT  data  to  be  on  the  order  of  25  to  40  km.  A  lower  resolution  could 
be  expected  from  the  GEOS-3  data  because  of  the  poor  noise  quality.  In 
addition  we  must  be  aware  of  the  variable  track  spacing  of  the  edtimeter  data 
which  can  vary  from  10  km  to  100  km  depending  on  geographic  area. 

In  addition  our  predicted  quantities  have  a  frequency  cut  off  implied  by 
the  data  selection  process.  For  most  of  our  predictions  we  selected  300  data 
points  in  a  4*x4*  block.  This  implies  one  data  point  per  26x26  km  block  which 
could  form  a  resolution  limit  of  our  data  which  would  be  on  the  order  of  0*.23 
which  would  correspond  to  a  spherical  harmonic  expansion  to  about  degree 
780. 


A  second  approach  to  this  problem  is  to  carry  out  a  power  spectrum 
analysis  of  the  data  using  Fast  Fourier  Transform  technique.  For  our  test 
eight  profiles  whose  length  is  60*  (i.e.  481  points)  were  chosen  in  different 
parts  of  the  oceans.  An  FFT  analysis  was  performed  by  T.  Engelis.  Adopting 
a  somewhat  arbitrary  criterion  that  99.7%  of  the  power  in  the  spectrum  defines 
the  maximum  wave  number  (m)  present  in  the  profile,  the  value  of  m  varied 
from  134  to  202  with  an  average  veilue  of  155  where  the  fundamental  frequency 
is  1  cycle/60  degrees.  Using  equation  cll  (in  Appendix  C)  these  wave  numbers 
correspond  to  spherical  harmonic  degrees  of  802,  1204,  and  928  respectively. 
In  terms  of  resolution  this  corresponds  to  24.9  km  to  16.5  kra,  or  an  average 
21.5  km  which  corresponds  to  0*.19.  (This  resolution  is  somewhat  smaller  than 
the  17  km  resolutions  from  Seasat  data  estimated  by  Brammer  and  Sailor 
(1980)).  On  the  basis  of  this  analysis  the  approximate  prediction  interval  would 
have  been  0*.19  so  that  using  0*.125  can  not  be  completely  justified.  However, 
it  is  clear  information  is  lost  if  a  0*.25  grid  is  used.  In  addition  the 
computation  of  the  values  on  the  0*.125  grid  has  been  done  in  a  consistent 
manner  considering  a  collocation  estimation  and  such  values  can  be  helpful  in 
areas  where  the  gravity  field  is  changing  rapidly. 


l*xl*  Mean  Anomaly  PeterminutionH  and  Comparisons  with  Terrestrial  Data 

Our  past  work  with  anomaly  recovery  from  altimeter  data  has  primarily 
been  for  l*xl*  anomalies  (Rapp  1979).  Many  of  the  computed  values  can  be 
compared  to  terrestrial  estimates  for  accuracy  checks  or  for  the  identification 
of  significant  errors  in  one  set  of  (primarily  terrestrial)  anomalies. 

Using  the  0*.125  gridded  data  a  set  of  l*xl*  anomalies  was  computed.  The 
computations  involved  the  straight  mean  of  all  avEulable  points  (usually  81)  that 


fell  in  and  on  the  borders  of  the  block.  The  disadvantage  of  this  procedure  is 
that  it  causes  a  small  error  correlation  between  adjacent  blocks  due  to  using 
some  of  the  same  data  points  in  the  mean  anomaly  estimations.  Between 
latitude  limits  *  72*,  37419  l*xl*  values  were  computed.  Of  these,  2364  were  on 
land  and  are  considered  unreliable.  The  more  reliable  means  were  obtained 
between  the  latitude  limits  of  *65*  where  there  were  34914  values  1961  of  which 
were  on  land.  The  accuracy  of  each  of  the  mean  l*xl*  anomalies  was  computed 
from  the  following  expression: 

n 

®i 

m(l*xl*)  =  -ii-! -  .  .  _i _  (22) 

n  n  9  97 


where  n  was  the  number  of  points  used  in  forming  the  mean,  mj  is  the 
standard  deviation  of  the  point  value,  and  2.27  was  an  empirical  factor  defined 
to  give  values  of  m(l*xl*)  that  were  consistent  with  more  rigorously  derived 
estimates.  No  mean  value  was  estimated  if  n  was  less  than  40.  To  identify 
such  blocks  a  m  value  of  999.0  was  assigned.  The  location  of  29880  l*xl* 
values  where  the  standard  deviation  is  <  6  mgals  is  shown  in  Figure  42. 


The  altimeter  derived  l*xl*  mean  anomalies  were  compared  to  our  must 
recent  set  of  l*xl*  anomalies  based  on  terrestrial  data  (the  January  1983  set) 
which  is  described  by  Rapp  (1983).  Several  comparisons  have  been  made  that 
vary  with  the  following  criteria:  maximum  standard  deviation  of  altimeter 
anomaly;  maximum  standard  deviation  of  terrestrial  anomaly  and  the  limit  on  the 
maximum  difference  allowed  in  the  comparisons.  The  results  of  these 
comparisons  are  given  in  Table  13.  No  comparisons  are  made  if  a  block  is  on 
land.  Comparisons  are  only  made  between  the  latitude  limits  *  65*. 


Figure  42 

Location  of  29880  Altimeter  Derived  l*xl*  Anomalies  with  Standard  Deviations  ^  6  mgals 


Table  13 

Comparison  of  Altimeter  Derived  I'xl*  Anomalies  to  Terrestrial  Estimates 
(January  1983  Data;  All  anomaly  values  are  in  mgals) 


Solution  Criteria 

SD  Alt  - 

_ 

12 

12 

12 

7 

SD  Terr  < 

— 

20 

15 

15 

15 

Diff  < 

— 

— 

— 

25 

25 

Quantity 

Num  of  Diff 

21933 

19419 

11463 

11183 

10139 

Min  Val  (Alt) 

-259 

-259 

-259 

-259 

Min  Val  (Terr) 

-282 

-282 

-282 

-282 

Max  Val  (Alt) 

167 

167 

167 

167 

Max  Val  (Terr) 

192 

192 

192 

192 

RMS  Val  (Alt) 

24.5 

24.5 

27.8 

27.5 

23.2 

RMS  Val  (Terr) 

27.2 

27.0 

29.7 

28.8 

24.5 

Mean  (A-T) 

0.1 

0.1 

0.5 

0.4 

0.4 

RMS  Diff 

*11.7 

*10.7 

*9.3 

*7.7 

*7.3 

RMS  SD  (Alt) 

5.4 

5.2 

5.4 

5.3 

4.8 

RMS  SD  (Terr) 

15.4 

14.2 

10.9 

10.9 

10.9 

Max  SD  (Alt) 

28 

12 

12 

12 

7 

Max  SD  (Terr) 

65 

20 

15 

15 

15 

Min  SD  (Alt) 

4 

4 

4 

4 

4 

Min  SD  (Terr) 

1 

1 

1 

1 

1 

Min  Diff  (A-T) 

-131 

-131 

-86 

-25 

Max  Diff  (A-T) 

130 

88 

88 

24 

Figure  43  shows  the  location  of  905  values  where  the  difference  between 
the  altimeter  derived  anomalies  (with  standard  deviation  less  than  or  equal  to 
12  mgals)  and  the  terrestrial  values  of  the  January  1983  tape  (with  no 
standard  deviation  limit  on  the  terrestrial  value)  exceeded  25  mgals.  A  total  of 
21054  anomalies  were  compared  in  this  analysis. 

Many  of  these  disagreements  are  due  to  errors  in  the  terrestrial  data.  We 
have  estimated  a  number  of  new  l*xl*  terrestrial  values  using  data  not 
avEiilable  in  the  January  1983  up  date.  A  summary  of  some  of  these 
comparisons  is  given  in  Table  14. 
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Figure  43 

Location  of  I'xl*  Anomalies  When  the  Differences  Between  the  Altimeter  Derived  Values  and  Terrestrial 

Estimates  Exceeds  25  mgals 


Table  14 

Selected  Comparisons  of  l*xl*  Anomalies  from  Altimeter  Estimates, 
January  1983  Data  Tape,  and  New  Terrestrial  Estimates 

(units  are  m^als) 


♦ 

X 

^6 

Altimeter 

January  83 

Difference 

New 

Estimate 

Source 

52 

175 

-22*9 

70*24 

-92 

-66*25 

A 

44 

350 

18*7 

-40*11 

58 

35*15 

B 

43 

350 

39*8 

-13*14 

52 

33*10 

B 

31 

292 

-17*4 

-105*13 

88 

-20*10 

C 

30 

292 

-18*4 

-73*15 

55 

-25*10 

C 

23 

268 

-8*5 

-75*15 

67 

-13*20 

D 

1 

6 

22*7 

108*15 

-86 

-5*10 

E 

-3 

327 

2*5 

58*16 

-56 

-10*10 

E 

-5 

147 

21*6 

-23*6 

44 

4*10 

F 

-6 

325 

0*5 

54*8 

-54 

-20*20 

E 

-15 

354 

6*5 

46*15 

-40 

-15*20 

E 

-31 

326 

16*6 

-55*24 

-71 

20*15 

E 

-32 

1 

7*4 

62*16 

-55 

22*10 

E 

-34 

316 

-17*5 

42*24 

-59 

-25*10 

E 

-38 

172 

9*4 

-49*14 

58 

11*14 

G 

-38 

173 

22*5 

-26*18 

48 

21*8 

G 

-38 

180 

29*5 

-11*21 

40 

35*8 

G 

-38 

181 

31*4 

-22*17 

53 

33*8 

G 

-38 

189 

-2*5 

-45*19 

43 

10*8 

G 

-41 

0 

-10*7 

65*20 

-75 

19*20 

E 

-41 

359 

-31*5 

72*25 

-103 

19*20 

E 

-51 

181 

-19*5 

64*19 

-83 

-18*7 

G 

The  sources  of  the  new  anomaly  data  are  as  follows: 

A;  Watts,  A.,  Gravity  Field  of  the  Northwest  Pacific  Ocean  Basin  and  It’s 
Margin,  Chart  MC-10,  Geological  Society  of  America,  Boulder,  Colorado,  1975 

B:  Torge,  W.,  et  al.,  6’xlO’  Free  Air  Anomalies  of  Europe  Including  Marine 
Areas,  Geophysical  Researchers,  Vol.  7,  Nos.  1/2,  1984 

C:  Watts,  A.,  (Ship  data,  private  communication,  1983) 

D:  Dehlinger,  P.,  Marine  Gravity,  Elsevier  Oceanography  Series,  1978 

E:  Kogan  M.G.,  and  A.  Watts,  Free  Air  Gravity  Field  of  the  South  Atlantic 
Ocean,  map  published  under  the  auspecies  of  the  USSR  Academy  of 
Sciences,  1983 
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F:  Watts,  A,  J.  Bodino,  C.  Bowin,  A  Geophysical  Atlas  of  the  East  and  Southeast 

Asian  Seas,  Free  Air  Gravity  Field,  Chart  MC-25,  The  Geological  Society  of 

America,  Boulder,  Colorado,  1978 

G:  Davey,  F.  and  A.  Watts,  Free-Air  Gravity  Filed  of  the  New  Zealand  Area, 

Chart  MC-48,  The  Geological  Society  of  America,  Boulder,  Colorado,  1983 

Inspection  of  Table  14  clearly  shows  that  the  new  altimeter  estimates  agree 
much  better  with  the  new  terrestrial  estimates  than  with  the  values  available  in 
the  January  1983  terrestrial  data  tape. 

If  we  consider  Table  13,  the  best  agreement  between  the  terrestrial  data 
and  the  altimeter  data  is  *7.3  mgals  which  is  for  the  case  when  unreliable 
values  in  both  data  sets  are  discarded.  This  value  is  substantially  smaller 
than  the  expected  value  of  *11.9  mgals  based  on  the  average  accuracy  of  the 
altimeter  anomalies  and  the  terrestrial  anomalies. 

0‘.5x0‘.5  Mean  Anomaly  Determination 

For  future  applications  in  the  development  of  high  degree  spherical 
harmonic  expansions  we  will  need  OtSxO’.S  mean  gravity  anomalies.  Such  values 
can  easily  be  computed  from  our  point  grid  data  in  a  manner  analogous  to  the 
l*xl*  mean  value  computation.  In  this  case  all  available  points  in  the  OtS  block 
were  meaned.  The  maximum  number  of  points  would  be  25.  The  minimum 
number  of  points  for  which  a  mean  was  determined  was  13.  The  0*.Sx0*.5  means 
were  determined  for  the  northern  hemisphere,  the  southern  hemisphere  and 
then  merged.  Some  basic  statistics  on  this  set  are  given  in  Table  15  for  two 
cases  selected  on  the  basis  of  the  anomaly  standard  devnation.  In  making  the 
statistics  we  have  no  easy  way  to  eliminate  land  values.  We  know  from  our 
l*xl*  data  set  that  there  were  2364  l*xl*  anomalies  on  land  or  about  9456 
0’.5x0*.5  mean  values.  The  RMS  value  of  *26.4  mgals  is  slightly  smaller  than  the 
RMS  value  of  the  point  anomalies  (*28.6  mgals)  and  higher  than  the  RMS  value 
of  all  the  l*xl*  mean  anomalies  (*22.8  mgals). 

The  0*.5  anomalies  computed  from  the  point  predictions  were  compared  to 
terrestrial  estimates  for  the  areas  described  earlier  in  this  report:  Japan 
Trench,  Puerto  Rico  Trench  and  the  Blake  Plateau.  (See  Tables  2,  3,  and  4). 
No  substantial  difference  was  found  between  the  new  results  and  the  ones 
described  in  these  tables. 

The  accuracy  estimates  of  the  0*.5  anomalies  were  computed  in  an  empirical 
way  similiar  to  that  done  for  the  l*xl*  anomalies.  The  overall  accuracy  of  *8 
mgals  may  be  somewhat  smaller  than  expected  •'"om  a  more  rigorous  treatment 
of  the  prediction  process. 


Table  15 

Basic  Statistics  on  0*.5x0*5  Mean  Anomalies  Derived  from  the  Altimeter  Data 

(units  are  mgals) 


Northern 

Hemisphere 

Southern 

Hemisphere 

Globally 

Nunbert 

59706 

89964 

149670 

RMS  Ag 

30.04 

23.6 

26.4 

Max  Ag 

319 

203 

319 

Min  Ag 

-286 

-241 

-286 

Max  SD 

55 

52 

55 

Min  SD 

4 

4 

4 

RMS  SD 

8.0 

8.0 

8.0 

Number! 

59346 

89481 

148827 

RMS  Ag 

29.8 

23.3 

26.1 

Max  Ag 

319 

193 

319 

Min  Ag 

-286 

-241 

-286 

Max  SD 

30.0 

30.0 

30.0 

Min  SD 

3.7 

4.4 

3.7 

RMS  SD 

7.5 

7.7 

7.6 

♦All  values  are  considered 

I  Only  values  where  the  standard  deviation  ^  30  mgals  are  considered 


In  making  these  statistics  we  have  no  easy  way  to  eliminate  land  values.  We 
know  from  our  l*xl*  data  set  that  there  were  2364  l*xl*  anomalies  on  land  or 
about  9456  0*.5x0*.5  mean  values.  The  RMS  value  of  »  26.4  mgals  is  slightly 
smaller  than  the  RMS  value  of  the  point  values  (*  28.6  mgals)  and  higher  than 
the  RMS  value  of  all  the  l*xl*  mean  anomalies  22.8  mgals). 

The  0*.5  anomalies  computed  from  the  point  predictions  were  compared  to 
terrestrial  estimates  for  the  areas  described  earlier  in  this  report:  Japan 
Trench,  Puerto  Rico  Trench  and  Blake  Plateau.  (See  Tables  2,  3,  and  4).  No 
substantial  differenace  was  found  between  these  new  results  and  the  ones 
described  in  these  tables. 

The  accuracy  estimates  of  the  0*.5  anomalies  was  computed  in  an  empirical 
way  similiar  to  that  done  for  the  l*xl*  anomalies.  The  overall  accuracy  of  *  8 
mgals  may  be  somewhat  smaller  than  expected  from  a  more  rigorous  treatment 
of  the  prediction  process. 


Potentiwl  fVwfficient  Behavior 


The  eighth  of  a  degree  data  set  can  provide  information  on  the  behavior  of 
potential  coefficients  used  to  describe  the  earth’s  gravitational  potential.  This 
is  done  by  examining  the  power  spectra  determined  in  local  areas  using  Fourier 
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analysis  procedures  (Forsberg,  1984).  Let  the  disturbing  potential  T  be 
represented  in  fully  nomalized  potential  coefficients  Cimi  Sin,  as  follows; 


V  • 


where 


V.  xi  -  I  IlmCsin?)  cos  mX;  m  ^  0 
’  ~  1  P4ni(sin?)  sin  mX;  n  <  0 


The  potential  degree  variances  of  T  are  designated  ^4! 


=  (?r  (fi'T.  ci. 


where  <n  is  formally  the  spectrum  at  radius  r.  The  units  of  at  will  be  in  some 
geopotential  unit  such  as  (1  kgal  m  =  Igpu).  The  unitless  degree  variances  are 
defined  as  ag: 


'•  -  c«. 


Comparing  (24)  and  (25)  we  have; 


/  fkMl* 

'4  =  Ir)  ‘ 


Gravity  anomalies  are  related  to  the  disturbing  potential  by  the  following 
spherical  approximation: 


Ag(r.  J,  X)  = 


Using  (23)  in  (27)  we  can  write: 


CiniYfni(^« 


(28) 


.8(r,  »,  X,  =^C(-0(S)'i: 

1  =  2  in=— < 


The  anomaly  degree  vairiancea  on  a  sphere  of  radius  a  would  be: 


C| 


(29) 


or 


C| 


(l-i)> 

i:'r  ‘  <"i 


(30) 


Since  kM/a*  =  r  (an  average  value  of  gravity  taken  as  979.8  gals)  we  can  also 
write  (29)  as: 


C|  =  7*(<-»)*V|  (31) 

We  can  also  talk  of  spatial  cov8u*iance  function  K(f,  rp,  rQ)  involving  the 
disturbing  potential.  In  the  planar  case  we  can  write  K(8(  zp,  zq)  where 
8  5  -fR  and  z  are  elevations  above  a  reference  plane.  Given  the  covariance 
function  K(8)  (we  assume  zp  =  zq  =  0)  the  power  spectrum  (or  power  spectral 
density)  can  be  written  as  (Porsberg,  eq  (2.28),  1984): 


♦^(«)  =  27t/*  sK(s)Jo(«s)ds  (32) 

where  u  is  a  continuous  frequency  function  and  Jo  is  the  Bessel  function  of 
zero  order.  The  integral  defines  the  Hankel  transform  (♦tt("))  of  K(s) 
(Papoulis,  p.  140,  1981).  This  is  equivalent  to  a  Fourier  Transform  when 
circular  symmetry  is  assumed,  as  we  do  for  covariance  functions.  If  K(s)  is 
given  in  potential  units  squared  (e.g.  gpu^)  and  if  s  is  given  in  angular 
degrees  the  units  of  ♦tt(")  be  gpu’deg*.  Forsberg  (ibid)  shows  that: 


=  (*  +  ’i)  ^  7^  ♦Tt("  =  (33) 


where  R  is  a  mean  earth  radius.  The  anomaly  degree  vsuriances,  ci,  can  also  be 
found  from  the  psd  of  the  anomaly  covariance  function: 


C| 


(I  +  X) 


(34) 
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ms 


where  the  units  of  ♦gg  would  be  mgal*deg*  if  the  anomalies  are  in  mgal  and  s 
is  in  degrees.  The  value  of  <ri  can  be  found  using  (29)  and  (34): 


~  (<-!)*  ♦««[  H  J 


(35) 


Assuming  ♦gg  is  given  in  (mgal  deg)^  and  1*  angular  degree  corresponds  to 
111  km,  (29)  Decomes: 


<T|  =  5.02x10  (♦gg  (-^)(mgal  deg)*] 


(36) 


Our  problem  now  is  to  estimate  the  anomaly  covariance  function  and  then 
compute  its  pad  and  then  compute  from  (36).  This  procedure  can  be  carried 
out  in  local  areas  using  our  eighth  degree  data  set.  We  have  done  this  using 
a  program  prepared  by  Forsberg  that  takes  an  input  grid  and  carries  out  a  2 
dimensional  Fourier  analysis.  The  individual  spectra  are  then  averaged  to 
obtain  a  specific  set  of  values  for  ♦tt  ♦gg. 

For  our  application  we  decided  to  compute  <r|  for  15*xl5*  areas  using  the 
eighth  degree  data.  To  reduce  spectral  leakage  we  removed  from  our  anomalies 
the  anomalies  implied  by  the  GEM9  potential  coefficients  taken  to  degree  20.  If 
NN  is  the  number  of  points  in  one  direction  of  the  grid  with  a  spacing  of  DPI 
(in  degrees)  the  frequency  spacing,  Af  is: 


Af  = 


(NN)(DFI) 


(37) 


The  frequencies  that  can  be  determined  will  be  iAf  where  1  <  i  (NN/2).  The 
corresponding  spherical  harmonic  degrees  are: 


li  =  i  360*  Af 


(38) 


Equation  (38)  is  the  same  as  (cll)  in  Appendix  C.  In  our  specific  case 
NN  =  481,  DPI  =  0M25  so  that  <  =  6i.  The  highest  spherical  harmonic  degree  of 
possible  interest  could  be  1440  (i.e.  6(481/2)).  However  earlier  discussion 
showed  that  we  should  restrict  our  analysis  to  maximum  <'s  on  the  order  of 
900. 


The  first  set  of  spectrum  calculations  were  made  for  the  eight  geographic 
areas  shown  in  Pigure  44  where  there  are  a  total  of  104  spectra. 
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Figure  44 

Location  of  Areas  in  Which  Local  Spectra  Were  (Tomputed.  Number  Below  Letter  is  the  Number  of  Spectra 

Computed  for  Bach  Area. 


The  spectra  from  each  15*xl5*  run  were  computed  and  the  results  menned. 
The  results  for  each  area  are  plotted  in  Figure  45.  Also  plotted  in  this  figure 
is  the  spectrum  due  to  Kaula’s  Rule  on  the  potential  coefficient  decay  and  the 
spectra  implied  by  the  anomaly  degree  variance  model  of  Tscherning  and  Rapp 
(1974).  The  Kaula  rule  spectrum  is: 


The  Tscherning/Rapp  spectrum  is: 

_  A _ (1-1)  I+, 

*  ~  (<-1)^7*  (l-2)(4+24)  ® 


(39) 


'40) 


where  A  =  425.28  mgaP  and  s  =  0.999617. 

The  spectrum  with  the  most  power  is  that  of  Area  G.  This  area  is 
dominated  by  anomalies  associated  with  trenches.  Area  H  is  the  area  that  has 
the  least  power.  This  area  is  characterized  by  small,  slowly  changing,  gravity 
anomalies.  The  other  areas  have  spectra  that  fall  between  the  two  extremes. 
The  dispersion  between  the  spectra  is  the  greatest  out  to  about  degree  360 
after  which  six  of  the  eight  areas  have  spectra  that  agree  quite  well.  The 
spectra  implied  by  the  two  models  represent  well  the  actual  spectra  in  area  B 
and  G  out  to  about  degree  300  after  which  the  models  contain  too  much  power. 

Next  the  104  spectra  were  divided  into  three  categories  that  were  selected 
to  describe  the  variability  of  the  anomaly  field.  To  do  this  the  range  of  the 
standard  deviations  (SD)  of  the  anomalies  with  respect  to  GEM9  was  divided  by 
three  to  yield  the  following  categories:  Rough  Area  (43  <  SD  <  59  mgals);  Mild 
Area  (23  <  SD  ^  43  mgals);  Smooth  Area  (4  *  SD  <  23  mgals).  The  number  of 
spectra  in  each  area  was:  Rough  (9),  Mild  (17),  Smooth  (78).  The  mean  SD  for 
each  area  was:  Rough  (*  52.1  mgals);  Mild  29.5  mgals);  Smooth  (*12.9 
mgals).  The  location  of  the  rough,  mild,  and  smooth  spectra  is  shown  in  Figure 
D1  in  Appendix  D. 

The  spectra  from  each  group  were  then  plotted  as  shown  in  Figure  46 
along  with  the  Kaula  and  Tscherning/Rapp  models.  The  mean  of  the  104 
spectra  is  also  shown.  It  is  clear  that  the  mean  spectrum  is  considerably 
smaller  than  expected  from  the  two  models.  For  more  meaningful  comparisons 
we  computed  the  RMS  coefficient  (7)  from  the  various  spectra.  We  have 


5f|(C,  S) 


(41) 


Values  of  are  compared  in  Table  16. 


Table  16.  Comparison  of  RMS  Coefficient  Variation 
(All  Values  xlO‘‘) 


Degree 

Mean  from  Spectra 

Kaula’s  Rule 

Ratio 

96 

50.6 

108.5 

2.15 

216 

14.6 

21.4 

1.47 

312 

7.5 

10.3 

1.37 

408 

4.3 

6.0 

1.41 

504 

2.6 

3.9 

1.49 

600 

1.7 

2.8 

1.65 

696 

1.0 

2.1 

1.98 

792 

0.7 

1.6 

2.35 

912 

0.4 

1.2 

2.85 

We  recognize  in  these  spectra  that  our  computed  values  are  only 
representative  of  the  ocean  regions  of  the  world.  Forsberg  (1984)  has 
computed  spectra  for  four  land  areas  in  the  United  States  using  gravity 
anomalies  given  on  a  2’x2’  grid.  The  four  areas  studied  by  Forsberg  were  in 
Ohioi  Colorado,  California  and  New  Mexico.  The  areas  were  classified  as  follows: 
Ohio:  lowland/moderately  hilly;  Colorado:  alpine  area  west  of  Denver  and 

Colorado  Springs;  California:  mixed  area  containing  the  central  parts  of  the 
Sierra  mountains  and  a  pau-t  of  the  California  Valley;  New  Mexico:  4*x4*  8u*ea 
covering  most  of  the  state.  The  anomalies  were  first  corrected  by  removing 
the  Bouguer  anomalies  implied  by  the  Dec  1981  potential  field  to  180  and  a 
spherical  harmonic  expansion  of  the  topography  to  degree  180.  A  mean  curve 
was  estimated  form  Figure  20  of  Forsberg  (1984)  and  plotted  in  Figure  46.  The 
curve  falls  near  the  mean  and  mild  curves  estimated  from  the  ocean  data. 

If  the  topographic  effects  had  not  been  removed  from  the  Colorado, 
California  and  New  Mexico  data,  the  spectra  would  be  more  varied  and  would 
show  more  power  than  the  spectra  computed  from  the  corrected  anomalies. 
This  is  clearly  seen  from  Figure  21  of  Forsberg  (ibid). 

Since  the  oceans  cover  70X  of  the  earth’s  surface,  and  since  rugged 
topographic  areas  are  not  the  normal  continental  structure,  the  mean  global 
spectrum  could  be  represented  by  a  spectrum  somewhat  between  the  mean  and 
mild  area  spectra  shown  in  Figure  46.  One  could  find  improved  spectrum 
models  by  fitting  to  these  curves. 

A  simple  computation  assumes  the  potential  power  spectrum  may  be 
represented  in  the  general  form: 


where,  for  example,  B  is  3  for  the  Kaula  rule.  Values  of  B  have  been 
determined  from  the  data  shown  in  Figure  45  with  selected  results  shown  in 
Table  17. 


Table  17.  Values  of  B  in  the  Potential  Decay  Law 


Data 

REinge  of  < 

B 

Forsberg 

200 

to 

500 

3.38 

600 

to 

1455 

3.74 

200 

to 

1455 

3.65 

Average  of 

24 

to 

500 

2.76 

104  Spectra 

24 

to 

792 

3.63 

24 

to 

960 

4.04 

100 

to 

700 

4.56 

100 

to 

792 

3.77 

500 

to 

792 

4.00 

792 

to 

1440 

5.69 

From  the  Forsberg  data  we  see  an  early  decay  of  3.38  which  increases  to 
3.74  for  degrees  600  and  above.  In  considering  the  altimeter  data  we  must 
recall  that  the  higher  frequencies  have  been  smoothed  by  the  collocation 
process  and  that  the  analysis  on  page  72  gave  the  highest  estimable  degree  to 
be  on  the  order  of  <=950.  The  slopes  shown  for  the  altimeter  analysis  in  Table 
17  indicate  the  more  rapid  decay  at  the  higher  frequencies.  This  behavior  can 
be  caused  by  several  factors:  1.  the  effect  of  the  smoothing  caused  by  the 
collocation  process  and  2.  the  fact  that  the  density  anomalies  are  farther  away 
from  the  surface  in  the  ocean  case  than  in  the  land  case.  Based  on  the 
comparison  with  the  slopes  of  the  Forsberg  data,  we  might  conclude  that  the 
altimeter  derived  spectra  shown  in  Figure  46  are  reasonable  up  to  degree  800 
after  which  they  probably  decay  too  rapidly. 


Anomaly  Degree  Variances 

An  alternative  to  the  analysis  of  potential  coefficient  behavior  is  the  study 
of  anomaly  degree  variances.  Such  values  can  be  computed  from  the  potential 
coefficient  spectrum  using  equation  (31).  The  anomaly  degree  variance  implied 
by  the  Kaula  Rule  is: 

ci  =  mgal*  (43) 


The  Tscherning/Rapp  model  was  given  in  terms  of  anomaly  degree  variances 


A0-A16C  591 
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DETAILED  GRAVITV  ANOMALIES  AND  SEA  SURFACE  HEIGHTS 
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DEPT  OF  GEODETIC  SCIENCE  AND  SURVEVI.  .  R  H  RAPP  AUG  85 
0SU/DGSS-TS5  HFGL-TR-85-ei91  F/G  8/5 


where  A,  B,  and  s  have  been  defined  after  equation  (40).  To  be  exact,  the  C| 
values  defined  by  (44)  refer  to  a  mean  earth  radius  of  6371  km  instead  of  an 
equatorial  radius.  We  will  ignore  the  distinction  here. 

Table  18  contains  the  anomaly  degree  variances  from  the  eight  geographic 
areas  discussed  earlier.  We  have  also  included  the  C|  values  computed  from 
the  two  models  given  in  (43)  and  (44).  And  we  also  have  computed  the  ci 
values  from  the  Rapp  (1981)  expansion  to  degree  180;  the  GEMIOC  potential 
coefficients  to  degree  180;  and  from  a  set  (Bl)of  potential  coefficients  to  degree 
250  estimated  by  Hajela  (1984)  in  an  "optimal"  way  using  the  same  adjusted 
anomalies  as  Rapp  (1981).  All  values  are  given  in  Table  18. 


Table  18  Anomaly  Degree  Variances  by  Geographic  Areas,  From  Models,  and 
Potential  Coefficient  Sets.  Units  are  mgal^. 


48 

72 

96 

Degree  ( t ) 
120  144 

168 

192 

216 

240 

Kaula 

3.96 

2.61 

1.97 

1.58 

1.32 

1.13 

0.99 

0.88 

0.80 

Tsch/Happ 

5.92 

4.37 

3.15 

2.84 

2.41 

2.09 

1.84 

1.64 

1.47 

Rapp(1981) 

2.64 

3.06 

2.45 

2.60 

2.13 

1.88 

— 

— 

— 

Hajela(1984) 

2.33 

2.78 

2.11 

2.17 

1.57 

1.08 

0.66 

0.48 

0.29 

CEMIOC 

2.62 

2.10 

1.56 

1.33 

1.10 

1.10 

— 

— 

— 

Region  A 

0.44 

0.49 

0.50 

0.61 

0.61 

0.63 

0.50 

0.56 

0.46 

B 

1.10 

1.11 

1.34 

1.61 

1.89 

1.56 

1.59 

1.41 

1.24 

C 

1.33 

0.98 

0.94 

0.68 

0.84 

0.66 

0.58 

0.48 

0.46 

0 

0.47 

0.50 

0.41 

0.47 

0.50 

0.53 

0.51 

0.44 

0.38 

E 

0.67 

0.45 

0.33 

0.34 

0.32 

0.43 

0.38 

0.35 

0.39 

F 

0.56 

0.40 

0.36 

0.31 

0.30 

0.29 

0.29 

0.30 

0.28 

G 

2.53 

2.60 

2.37 

2.30 

1.98 

2.29 

2.01 

1.65 

1.62 

H 

0.17 

0.11 

0.09 

0.08 

0.08 

0.09 

0.12 

0.12 

0.11 

Rough 

7.33 

11.32 

12.26 

11.27 

12.70 

11.10 

8.06 

6.24 

4.53 

Mild 

1.51 

2.02 

2.43 

2.94 

2.75 

2.74 

2.49 

2.08 

1.79 

Smooth 

0.36 

0.25 

0.20 

0.19 

0.19 

0.21 

0.21 

0.21 

0.21 

Mean  of  104 

0.59 

0.48 

0.43 

0.42 

0.43 

0.44 

0.44 

0.41 

0.38 

In  comparing  the  results  we  must  keep  in  mind  that  the  two  models  and 
the  spherical  harmonic  sets  represent  the  earth  as  a  whole  while  the  other 
values  eu*e  representative  of  the  ocean  regions  alone.  At  degree  72,  the  value 
of  c|  should  be  between  2  to  3  mgal’  as  inferred  from  the  spherical  harmonic 
expansions.  However,  from  the  ocean  spectra  ci  is  on  the  order  of  0.5  mgal^. 
In  principle  the  differences  might  be  caused  by  the  land  data  which  is 
included  in  the  spherical  harmonic  expansions.  Whether  a  factor  of  5  (at 
degree  72)  between  the  two  types  of  solutions  is  reasonable  is  not  clear. 

Estimates  of  ci  can  be  obtained  from  Figures  20  and  21  of  Forsberg  (1984). 
At  degree  200  C|  is  approximately  0.5  mgal^  when  using  terrain  corrected  data 
and  1.9  mgaP  when  the  non-corrected  data  is  used.  The  0.5  mgal^  is 


consistent  with  the  value  of  0.4  mgaP  found  from  the  mean  of  our  104  profiles. 
The  increase  of  ci(200)  to  1.9  mgfal^  indicates  the  sensitivity  of  this  problem  to 
the  terrain. 


The  major  question  from  this  discussion  relates  to  the  substantial 
difference  between  the  C|  values  derived  from  our  local  spectra  and  the  global 
spherical  harmonic  implied  values.  Is  the  difference  due  to  the  topography? 
Or  is  it  due  to  some  other  factor  in  the  analylsis? 


Geoid  Undulation  Information  by  Degree  Range 

The  geoid  undulation  is  related  to  the  disturbing  potential  as  follows: 


N  =  T/7 


(44) 


where  7  is  a  value  of  gravity. 

Using  (23)  (with  r  =  R)  the  geoid  undulation  at  degree  <  is: 

I 

N|(R,  ^  (I)  'L  C4ai^4iii(^»  (45) 

oi=— I 

The  unculation  variance,  per  degree,  is: 

..<N,  =  (gni)“ .. 

If  we  let  7  =  kM/R»  and  R  =  a  we  have: 

<"l(N)  =  R*v«  (47) 

The  geoid  undulation  variance  between  degrees  li  and  1 2  would  be: 


AN»(I,, 


(48) 


The  value  of  <it  has  been  estimated  from  our  altimeter  derived  anomaly  data 
using  equation  (36).  It  can  be  computed  from  the  models  given  by  equation 
(39)  and  (40),  or  alternate  models  such  as  given  by  Wagner  and  Colombo  (1979). 
Brammer  and  Sailor  (1982)  have  analyzed  Seasat  altimeter  data  and  terrestrial 


gravity  data  from  ship  measurements  to  produce  a  power  law  decay  model  for 
the  potential.  This  model  has  been  used  to  compute  the  geoid  undulation 
information  in  several  degree  ranges.  These  values  are  shown  in  Table  18 
along  with  values  from  our  rough,  mild,  smooth  areas  and  from  our  average 
spectrum,  as  well  as  from  the  Kaula  rule,  the  Tscherning/Rapp  model,  and  the 
Wagner/Colombo  model. 

We  can  see  the  variability  of  the  results  depending  on  the  area 

characteristic  (rough,  mild,  or  smooth).  The  total  average  results  are  close  to 
the  results  for  the  smooth  areas,  because  the  latter  dominates  in  the  ocean. 
The  results  from  the  Tscherning/Rapp,  Wagner/Colombo,  and  Kaula  models  have 
more  power  than  seen  from  the  local  spectra  information.  (The 
Tscherning/Rapp  and  the  Wagner/Colombo  models  are  very  similiar  above 

degree  100).  The  Brammer/  Sailor  results  agree  well  with  our  results  (i.e.  total 

average)  in  the  ranges  181-360,  361-720,  and  721-1000.  Because  of  the 

geographic  dependence  on  the  arc  selection  or  ship  tracks  the  results  of 
Brammer/Sailor  will  depend  on  the  location  of  the  profiles  that  were  analyzed. 

Table  Gooid  Undulation  Information  Between  Selected  Spherical  Harmonic 

Degrees  (Units  are  meters.) 


Range  of 

1 

Rough 

Area 

Mild 

Area 

Smooth 

Area 

Total 

Average 

Tscheming 

Rapp 

Wagner 

Colombo 

Kaula 

Brammer 

Sailor 

50-100 

2.16 

.93 

.34 

.47 

1.44 

1.17 

1.12 

1.12 

101-180 

1.49 

.73 

.19 

.28 

0.71 

0.68 

0.53 

0.49 

181-360 

0.76 

.45 

.16 

.21 

0.42 

0.40 

0.31 

0.26 

361-720 

0.20 

.14 

.09 

.10 

0.20 

0.19 

0.15 

0.12 

721-1000 

0.04 

.03 

.02 

.03 

0.08 

0.07 

0.06 

0.04 

Maximum  Anomaly  and  Sea  Surface  Height  Gradients 

We  now  Mrish  to  identify  those  areas  that  have  the  largest  changes  in  the 
anomalies  and  sea  surface  heights.  We  first  computed  the  east/west  and 
north/south  gradients  using  the  0*.5x0*.5  mean  anomalies  or  sea  surface  heights. 
The  total  gradient  was  computed  by  combining  the  individual  components.  The 
maximum  values  were  identified  and  then  maps  of  each  area  were  made  and 
specific  gradients  were  computed  based  on  the  point  data  interpolated  from  the 
contour  maps. 

The  largest  anomaly  gradient  was  found  over  the  Kuril  Trench  in  the 
northwest  Pacific  Ocean.  This  gradient  was  6.9  mgal/km  near  ♦  =  47*.5, 
X  =  154*5.  The  anomaly  and  sea  surface  height  maps  in  this  area  are  shown  in 
Figure  47.  There  are  a  number  of  areas  where  the  anomaly  gradient  reaches  5 
mgal/km  or  somewhat  higher.  Such  gradients  should  be  regarded  as  an 
average  gradients  because  of  the  smoothing  that  has  taken  place  in  the  data 
analysis. 
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The  maximum  sea  surface  height  gradient  was  0.28  m/km.  The  location  of 
this  gradient  was  near  ♦  =  44t5,  X  :  148t5.  This  is  also  across  the  Kuril  Trench 
but  not  as  same  location  as  where  the  anomaly  gradient  is  a  maximum.  The 
above  gradient  corresponds  to  a  deflection  of  the  vertical  of  58". 

It  would  be  of  interest  to  prepare  a  map  showing  these  gradients  in  some 
global  fashion.  However,  the  only  global  set  of  gradients  we  have  is  based  on 
0*5  mean  values  which  may  not  show  enough  detail.  This  needs  to  be  tried. 
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Point  Anomaly  Variance 


We  saw  that  the  root  mean  square  value  of  the  predicted  point  anomalies 
(for  values  with  the  standard  deviation  less  than  or  equal  to  30  mgals)  was 
*  28.6  mgals  or  a  variance  of  818  mgal^.  This  value  does  not  reflect 
contributions  from  spherical  harmonic  degrees  above  degree  900 
(approximately).  An  estimate  of  this  contribution  can  be  made  using  anomaly 
degree  variance  models.  A  very  approximate  estimate  (based  on  a  modified 
Tscherning/Rapp  model)  is  200  mgal^.  Thus  we  have  a  point  anomaly  variance 
in  the  ocean  areas  of  about  1000  mgal’  or  *  32  mgal.  This  is  considerably 
smaller  than  the  1795  mgals^  (*  42  mgal)  of  the  original  Tscherning/Rapp 
model. 

For  other  comparisons  of  the  point  variance  we  used  the  results  of  Schwarz 
and  Lachapelle  (1980)  for  Canada  where  they  found  an  average  variance  of  1115 
mgaP  (*  33  mgals)  which  agrees  well  with  the  altimeter  result. 

We  also  calculated  the  variance  from  the  6'xlO’  anomalies  described  by 
Torge  et  al  (1983)  for  most  of  western  Europe.  In  this  case  we  have  a  variance 
of  1218  mgal*  (*  35  mgals).  Again  this  value  is  reasonably  consistent  with  the 
altimeter  value  and  the  Canadian  result. 

It  seems  clear  from  this  brief  analysis  that  the  global  anomaly  variance  of 
1795  mgal*  computed  by  Tscherning  and  Rapp  in  1974  is  too  high  with  a  more 
reasonable  value  being  around  1100  mgal*. 


Summary 

In  1982  we  completed  an  analysis  of  Seasat  altimeter  data.  This  analysis 
included  the  adjustment  of  the  Seasat  arcs  to  remove  orbit  error  using  a 
crossing  arc  technique.  This  data  set  was  used  to  estimate  l*xl*  mean  gravity 
anomalies  and  sea  surface  heights  (approximately  geoid  undulations).  In  1983 
Liang  carried  out  an  adjustment  that  merged  the  Geos-3  altimeter  data  to  the 
Seasat  data.  This  combined  data  set  contained  approximately  6  million  altimeter 
measurements.  This  data  set  was  analyzed  in  this  paper  for  an  improved 
determination  of  the  gravity  field  in  ocean  area. 

The  first  step  taken  was  to  carry  out  tests  for  the  recovery  of  mean  and 
point  gravity  anomalies  using  the  combined  Geos-3/Seasat  data.  These  tests 
were  made  with  0*.5x0’.5  mean  anomaly  and  point  gravity  anomaly  predictions 
with  the  comparison  of  the  predicted  values  with  "ground  truth"  defined  by 
ship  gravity  measurements.  These  tests  examined  the  role  of  various  data 
selection  procedures  and  the  use  of  different  covariance  models  in  the  least 
squeures  collocation  prediction  method  used  to  recover  the  gravity  anomalies. 
These  tests  clearly  indicated  the  ability  of  the  altimeter  data  to  yield  reliable 
{*  8  mgals)  0*5  anomalies  and  promising  (^  10  mgEils)  point  anomcdies.  (Actual 
accuracy  estimates  will  depend  on  the  anomaly  behavior  in  an  area  and  the 
density  of  the  altimeter  data).  We  found  that  the  Tscherning/Rapp  covariance 
model,  with  appropriate  variance  scaling  to  the  altimeter  data  in  a  prediction 
region,  gave  somewhat  better  predictions  than  when  a  covariance  model  (with 
two  adjusted  parameters)  suggested  by  Jordan  (1972)  was  used.  In  both  cases 
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the  reference  field  was  defined  by  the  spherical  harmonic  expansion  to  degree 
180  of  Rapp  (1981).  It  became  clear  with  these  tests  that  the  primary  quantity 
to  be  predicted  should  be  values  of  the  anomalies  and  sea  surface  heights  on  a 
point  grid.  Although  such  values  appear  to  be  point  values,  they  are  not 
because  of  frequencies  missing  in  the  predicted  quantity  caused  by  the  use  of 
finite  data  having  noise. 

Many  tests  were  then  carried  out  to  develop  an  efficient  process  for  the 
production  computation  of  the  point  values  from  the  combined  altimeter  data 
set.  Decisions  were  made  to  predict  the  point  values  within  a  3*x3*  block  on  a 
0*.125x0*.125  grid  using  a  maximum  of  about  300  data  points.  Consequently  one 
data  selection  and  one  matrix  inversion  was  needed  for  each  3*x3*  area  in  the 
ocean.  In  order  to  save  computer  time  we  deleted  the  ability  to  compute  the 
rigorous  accuracy  of  the  predicted  quantity.  Such  a  decision  reduced  our 
computer  requirements  by  a  factor  of  8.  In  addition  our  data  selection  was 
made  from  a  thinned  data  base  created  by  selecting  a  single  data  point  for 
each  lO’xlO'  geographic  cell.  This  procedure  reduced  our  original  5.9  million 
point  data  base  to  a  1.1  million  point  set.  The  thinned  data  set  reduced  our 
data  selection  time  and  reduced  track  error  effects  since  points  were  not 
systematically  selected  from  a  single  track.  However  this  thinned  data  set,  and 
the  decision  to  use  only  300  data  points  per  block  prediction,  does  reduce  the 
resolution  of  our  prediction  over  the  case  of  edl  data  being  considered  with  a 
large  number  of  data  points  being  selected  for  the  prediction.  The  latter 
procedure  causes  other  problems,  notably  incrsased  data  selection  times  and 
matrix  instabilities. 

Before  the  production  predictions  were  started  a  number  of  Seasat  arcs 
were  re-adjusted  to  reduce  crossover  discrepancies.  This  was  done  by 
breaking  the  arcs  into  two  components  and  adjusting  each  arc  separately.  In 
our  actual  computations  we  still  saw  examples  of  remaining  track  error  which 
easily  showed  up  in  anomaly  contour  maps.  In  most  of  these  cases  the  bad  arc 
was  identified  and  the  prediction  repeated  with  the  bad  arc  removed. 


We  also  considered  the  role  of  sea  surface  topography  (SST)  in  our 
predictions.  Our  initial  plan  was  to  remove  SST  for  our  altimeter  observations 
using  a  specific  SST  model.  We  decided  not  to  do  this  because  such  models 
may  not  be  representative  of  the  ocean  surface  at  the  time  of  the  altimeter 
measurement  (or  general  time  of  the  altimeter  measurements),  and  because  we 
did  not  want  our  results  to  be  dependent  on  a  specific  model.  Therefore  the 
predicted  sea  surface  heights  refer  to  a  mean  sea  surface  defined  by  the 
average  time  of  the  altimeter  measurements  which  were  in  the  1975-1978  time 
period.  The  values  contain  sea  surface  topography  effects.  The  gravity 
anomalies  that  were  predicted  will  contain  errors  caused  by  not  removing  SST 
from  the  altimeter  measurement.  On  a  point  basis,  this  error  was  found  to  be 
*  2  mgal.  This  error  is  also  to  be  considered  long  wavelength  reflecting  the 
structure  of  long  term  SST. 


A  total  of  2.4  million  point  values  were  predicted.  Of  these  0.2  million  had 
standard  deviations  greater  than  *  30  mgals.  (The  standard  deviations  were 
estimated  considering  data  density  and  the  roughness  of  the  anomaly  field; 
they  were  not  rigorously  computed.)  The  meucimum  anomaly  was  396  mgal,  the 
minimum  value  was  -300  mgals,  and  the  root  mean  square  point  value  was  *  28.6 
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royals.  The  eijChth  of  a  degree  data  set  was  used  to  produce  a  number  of 
anomaly  and  sea  surface  height  (SST)  maps  that  showed  the  high  frequency 
information  in  the  predicted  quantities,  as  well  as  the  consistency  of  the 
results.  The  maps  clearly  show  the  trenches,  seamounts,  fracture  zones 
associated  with  major  features  of  the  ocean  floor.  A  number  of  seamount 
signatures  are  clearly  seen  for  which  no  bathymetric  data  apparently  exists. 
Contour  maps  revealed  very  large  anomaly  and  SST  gradients.  The  largest 
anomaly  gradient  was  found  to  be  6.9  mgal/km  over  the  Kuril  Trench  in  the 
northwest  Pacific  Ocean.  In  this  area  we  found  an  anomedy  change  of  345 
mgals  in  just  50  km. 

The  frequency  content  of  the  point  predictions  was  examined  by  carrying 
out  a  power  spectrum  analysis  of  eight  long  profiles.  This  analysis  indicated 
that  the  average  resolution  of  the  predicted  data  was  0tl9  or  21  km.  This  is 
somewhat  lower  resolution  than  implied  by  the  prediction  grid  interval  of 
0*125. 

The  point  anomaly  values  were  used  to  form  l*xl*  and  0*.5x0*.5  mean 
anomalies  by  the  averaging  of  the  point  values  in  the  specific  block.  Accuracy 
estimates  for  each  block  were  computed  in  an  approximate  way.  A  total  of 
37419  l*xl*  values  were  estimated.  The  root  mean  square  standard  deviation 
was  *  5  mgals.  These  were  29880  values  having  a  standard  deviation  less  than 
or  equal  to  *  6  mgals.  The  l*xl*  data  were  compared  to  our  latest  terrestrial 
data.  Using  the  more  reliable  data  in  each  data  set,  the  RMS  difference  was  * 
11  mgals  for  19419  values.  If  one  deletes  residuals  greater  than  25  mgals  this 
difference  drops  to  *  7.7  mgals  for  11183  computed  values.  Comparisons  of 
some  anomaly  estimates  with  new  terrestrial  estimates  indicated  the  reliability 
of  the  altimeter  derived  values. 

There  were  148827  0*.5x0*.5  anomalies  estimated  that  had  a  standard 
deviation  of  *  30  mgals  or  smaller.  The  root  mean  square  anomaly  value  was 
*  26  mgals  while  the  root  mean  square  standard  deviation  was  *  8  mgals. 

Our  final  computation  with  the  data  was  to  calculate  the  potential 
coefficient  spectrum  implied  by  the  anomalies  for  104  15*xl5*  regions.  This 
was  done  using  a  procedure  developed  by  Forsberg  (1984).  These  spectra 
were  classified  in  terms  of  rough,  mild,  and  smooth  areas  through  the  anomaly 
variance  in  the  15*  block.  The  smooth  areas  clearly  dominate  in  the  ocean 
areas  (78  smooth,  9  rough,  17  mild).  Vfe  found  that  the  average  spectrum  was 
substantially  smaller  than  implied-  by  Kaula  Rule  for  potential  coefficient  decay 
or  implied  by  the  Tscherning/Rapp  anomaly  degree  variance  model.  However 
the  average  curve  agreed  well  with  the  curve  from  Forsberg  based  on 
topographically  reduced  land  data. 

The  point  anomaly  variance  implied  by  the  predicted  altimeter  anomalies  was 
about  1000  mgal^  after  considering  the  influence  of  neglected  high  frequency 
terms.  This  is  considerably  less  than  the  1795  mgal'  point  variance  computed 
by  Tscherning/Rapp  in  1974.  However  we  saw  that  this  smaller  variance  is 
more  consistent  with  the  anomaly  variance  found  in  Canada  and  in  Europe.  We 
also  found  that  the  anomaly  degree  variances  implied  by  our  spectrum 
computation  were  considerably  smaller  at  around  degree  168  (say)  than  found 
from  spherical  harmonic  expansions  out  to  degrees  180  or  240.  One  would 
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expect  these  expansions  to  yield  higher  values  due  to  their  use  of  land  data  in 
their  solution.  However  it  seems  unreasonable  that  the  land  data  could  cause 
the  values  to  differ  by  a  factor  of  four  or  more. 

The  results  of  this  study  provide  a  detailed  gravity  anomaly  field  in  the 
ocean  eureas  between  latitude  limits  ^  65*  excluding  certain  small  sea  areas. 
This  data  can  be  a  basis  for  geophysical  or  geodetic  studies  that  require 
knowledge  of  the  detailed  gravity  field.  However  the  detail  that  we  have 
obtained  is  not  the  most  that  can  be  extracted  from  the  altimeter  data.  With 
additioned  computer  resources  more  data  points  could  be  used  in  the  prediction, 
and  more  detail  in  the  actual  field  could  be  extracted,  as  well  as  actual 
accuracy  predictions  made. 

A  fundamental  unresolved  question  in  understanding  our  results  relates  to 
the  accuracy  of  our  adjusted  altimeter  data.  Our  original  Seasat  adjustment 
fixed  one  long  arc  in  the  Atlantic  Ocean.  The  coordinate  frame  in  which  this 
(and  other  arcs)  are  defined  is  that  defined  for  orbit  determination  with  the 
release  by  JPL  of  the  geophysical  data  seconds  (GDR).  The  gravity  field 
(PGS-S3)  and  station  coordinates  associated  with  the  GDR  are  discussed  by 
Lerch  et  al,  (1982).  All  our  adjusted  values  refer  to  this  arc.  If  this  arc  is  in 
error  so  will  all  the  other  arcs  and  such  errors  will  propagate  into  the 
recovered  anomalies  and  sea  surface  heights.  Such  problems  will  be 
encountered  if  the  bias  and  tilt  adjustments  have  not  removed  most  of  the 
orbit  error.  From  the  comparisons  that  have  been  done  with  ship  data  we  have 
not  found  any  systematic  pattern  that  would  indicate  systematic  errors  caused 
by  poor  orbit  adjustment.  However  such  comparisons  may  not  reveal  errors  at 
the  several  mgal  or  several  decimeter  level.  Additional  study  is  needed  in 
orbit  adjustment  of  the  altimeter  data  and  the  effect  on  gravity  anomaly  and 
sea  surface  height  computation. 

We  also  need  to  study  the  smoothing  effects  introduced  in  our  predictions 
through  the  use  of  the  least  squares  collocation  procedure.  Do  the  smoother 
specira  and  the  smaller  anomaly  variances  reflect  reality,  or  is  there  an 
unaccounted  for  smoothing  that  is  caused  by  the  collocation  process? 

We  must  also  develop  better  models  for  sea  surface  topography  that  are 
valid  for  the  average  time  of  the  altimeter  measurements.  Our  recovered 
anomalies  are  contaminated  by  SST  effects  that  have  not  been  removed  from 
the  altimeter  measurements.  Such  errors  will  primarily  be  long  wavelength 
because  we  have  averaged  the  altimeter  measurements  over  several  years.  Can 
we  develop  a  procedure  to  correct  the  predicted  anomalies  or  must  we  first 
correct  each  altimeter  measurement? 

The  gravity  field  information  obtcuned  in  this  report  (and  studies  such  as 
carried  out  by  Sandwell  (1984),  Haxby  et  al.  1983  and  others)  clearly  reveal  the 
unique  role  that  satellite  altimetry  can  play  in  the  determination  of  the  gravity 
field  of  the  oceans.  However  as  our  accuracy  requirements  increase,  we  will 
find  it  necessary  to  seek  improvements  in  our  processing  algorithms,  or 
understanding  of  altimeter  satellite  orbits,  sea  surface  topography,  and  other 
such  effects  for  which  only  approximate  models  are  now  aveulable. 
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The  Seasat  adjustment  was  originally  carried  out  in  two  stepsi  a  primary 
adjustment  and  a  local  adjustment.  The  primary  adjustment  was  a  crossing  arc 
adjustment  of  long  arcs  to  obtain  a  global  network  to  which  the  remaining 
(and  mostly)  shorter  arcs  could  be  adjusted  in  five  geographic  regions.  The 
549  adjusted  primary  arcs  are  shown  in  Figure  Al.  Although  the  RMS 
crossover  discrepancy  after  the  adjustment  of  the  primary  arcs  was  t  28  cm 
(Rowlands  1981,  p.51)  there  were  53  arcs  where  the  crossover  discrepancies 
(along  the  arc)  were  greater  than  40  cm.  The  location  of  these  arcs  is  shown 
in  Figure  A2.  (Of  these  arcs  the  largest  RMS  crossover  discrepency  was  t61 
cm.)  The  reason  for  the  larger  RMS  is  primarily  due  to  the  length  of  the  arc 
being  too  long  for  the  radial  orbit  error  to  be  modeled  by  the  simple  bias  and 
tilt  term.  To  reduce  this  error  we  decided  to  readjust  these  arcs  by  splitting 
them  into  smaller  arcs. 

The  readjustments  were  regional  adjustments.  By  giving  the  size  uf  a 
region,  the  observations  of  these  53  arcs  were  collected  only  for  those 
observations  that  fell  into  the  region.  Such  observations  from  an  arc  were 
collected  to  form  a  cut  arc.  For  one  cut  arc,  the  crossovers  of  the  cut  arc  and 
the  remaining  primary  arcs  (excluding  the  53  arcs)  were  found  and  used  to 
readjust  the  cut  arc  by  a  crossover  adjustment.  In  other  words,  the 
readjustment  was  a  crossover  adjustment  and  was  separately  performed  in 
principle  for  each  cut  arc.  Practically,  all  cut  arcs  within  one  region  were 
adjusted  with  one  regional  adjustment  although  they  were  independent  of  one 
another. 

From  Figure  A2,  we  see  that  most  of  the  high  crossover  discrepancy  arcs 
were  located  in  the  Pacific  area,  therefore,  two  regional  adjustments  for  the 
Pacific  area  and  another  two  for  the  remaining  areas  were  performed.  The 
results  are  shown  in  Table  Al. 
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Figure  A1 

Distribution  of  the  SEASAT  Primary  Arcs. 
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Figure  A2 

53  Arcs  Having  Crossover  Discrepancies  Larger  Than  40 


Table  A1 

Locations  and  Statistics  of  Selected  Primary  Arc  Adjustment 


Region 

N.  Pacific 

S.  Pacific 

North 

South 

Total 

Latitude 

-2*<4<72* 

-72*i^*2* 

-2*^4i72* 

-72*^^<2* 

Longitude 

158*E<X-698*W 

158*E^X^98*W 

102‘W<X^162*E 

102*W<X<162*E 

no.  of 
cut  arcs 

23 

30 

16 

29 

98 

no.  of 

crossovers 

256 

523 

93 

362 

1234 

♦priori  RMS 
crossover 
discrepancies 
cm 

40 

36 

36 

31 

35.5 

♦♦posteriori 
RMS  cross¬ 
over  dis¬ 
crepancies 

cm 

22 

21 

23 

24 

22.3 

♦  after  the  primary  adjustment 
♦♦  after  the  readjustment 


After  the  readjustments,  the  data  sets  containing  the  primary  adjusted 
arcs  were  updated.  Only  the  adjusted  sea  surface  heights  of  the  readjusted 
arcs  were  replaced.  The  remaining  data  points  were  kept  the  same  as  before. 
With  the  updated  data  sets,  (J;25x0*.25  grid  collocation  predictions  were  carried 
out  and  the  results  compared  to  values  predicted  from  the  previously  adjusted 
Seasat  data.  The  irregular  track  patterns  that  were  present  in  the  latter 
results  were  no  longer  present  in  the  revised  data  set.  It  was  cleeu*  that  a 
substantial  improvement  had  been  made  with  the  revised  Seasat  data.  We 
should  note  that  no  readjustment  of  the  local  arcs  to  the  newly  adjusted 
Seasat  arcs  was  made.  In  principle  this  should  be  done  but  the  effect  in  not 
doing  this  is  expected  to  be  small. 


In  order  to  determine  the  geoid  from  altimeter  measurements,  an  estimate  of 
the  stationary  sea  surface  topography  (SST)  is  required.  The  stationary  SST 
represents  the  deviation  of  the  mean  sea  surface  from  the  geoid 
(appropriately  defined),  caused  mainly  by  oceanic  circulation  phenomena.  The 
geoid  in  turn,  has  to  be  defined  in  such  a  way  that  the  average  non  tidal  SST 
is  zero,  for  a  specific  epoch  of  measurement,  as  sampled  globally  in  the  oceanic 
regions  (Rizos,  1980). 

To  date,  only  estimates  of  the  large  scale  SST  are  available.  These 
estimates  have  been  computed  either  by  geodetic  techniques  (Engelis  1983, 
1985,  Engelis  et  al  1984,  Tai  and  Wunsch  1984)  or  by  traditional  oceanographic 
methods  (Levitus,  1982). 

Geodetic  techniques  require  an  independent  estimate  of  the  marine  geoid, 
which  is  provided  only  by  satellite  derived  gravity  fields.  These  gravity  fields 
have  accuracies  which  are  acceptable,  for  SST  determinations,  only  for  very 
long  wavelengths.  Consequently,  only  the  large  scale  SST  can  be  adequately 
determined  and  used  to  correct  the  mean  sea  surface  to  obtain  an  estimate  of 
the  geoid. 

In  recent  studies  by  Engelis  (1985)  long  wavelength  SST  was  determined  .i.s 
follows:  a  set  of  1*  gridded  sea  surface  heights  was  used  to  represent  the 

mean  sea  surface  as  realized  by  SEASAT.  The  geoid  was  implied  by  using  the 
GEML2  gravity  model  up  to  degree  6  (providing  an  accuracy  of  18  cm  in  geoidal 
heights).  Then  the  difference  of  the  two,  up  to  degree  6,  would  in  principle 
provide  the  large  scale  SST.  In  order  to  reduce  the  leakage  of  energy 
between  different  wavelengths  due  to  the  inconsistent  representation  of  the 
geoid  and  the  mean  sea  surface,  the  gravity  model  by  Rapp  (1981)  was  used 
from  degree  7  to  180,  to  create  geoidal  heights  with  wavelengths  comparable  to 
the  ones  of  the  sea  surface  heights.  Then,  raw  estimates  of  SST  were  obtained 
by  differencing  the  two  sets.  In  order  to  conform  with  the  definition  of  the 
geoid,  the  mean  value  of  these  estimates  was  obtained  and  a  new  data  set  was 
created  that  had  a  zero  mean.  In  this  set  zeros  have  been  used  for  land  and 
the  polar  regions.  This  global  set  has  been  harmonically  analysed  and 
harmonic  coefficients  up  to  degree  6,  representing  SST  effects,  were  retained. 
Synthesis  of  these  coefficients  provided  maps  of  SST  up  to  degree  6  that  can 
be  found  in  Engelis  (1985).  Errors  that  are  to  be  associated  with  such  a 
determination  of  the  SST  include  errors  in  the  GEML2  model  up  to  degree  6, 
very  long  wavelength  orbital  errors  remaining  from  the  crossover  adjustment 
and  leakage  of  energy  between  coefficients  due  to  errors  in  the  Rapp  model 
and  shorter  wavelength  orbital  errors.  The  estimated  accuracy  of  the  SST 
estimates  is  on  the  order  of  20  to  25  cm. 

Sea  surface  topography  has  also  been  determined  by  oceanographic 
methods.  These  methods,  being  unable  to  use  an  estimate  of  the  geoid,  define 
a  deep  surface  as  a  level  of  no  motion  (equipotential  surface)  and  then  they 
determine  the  mean  annual  anomaly  of  geopotential  thickness  of  the  layer 
between  that  surface  and  the  surface  of  the  oceans,  using  geostrophy 
assumptions.  The  data  used  are  salinity,  temperature,  and  dissolved  oxygen 
observations.  The  most  recent  determination  of  such  a  dymanic  sea  surface 
topography  was  made  by  Levitus  (1982).  A  2000  dbar  deep  surface  was  used 
as  a  reference.  The  data  that  Levitus  has  used  were  obtained  by  the  National 


Oceanographic  Data  Center.  This  data  file  represents  adl  observations  of 
salinity,  temperature,  and  dissolved  oxygen  taken  over  several  decades.  A 
first  analysis  made  by  Levitua  indicated  severed  problems.  The  most  important 
were  regional  biases  of  data  and  lack  of  data  in  extended  regions.  Moreover 
there  were  representative  problems  of  the  data,  since  observations  were  not 
synoptic  but  scattered  with  respect  to  time  (with  the  exception  of  a  few  limited 
areas)  and  so,  the  results  cannot  in  a  strict  sense  be  considered  a  true  long 
term  average. 

After  the  initial  stability  and  statistical  checks  to  eliminate  spurious 
observations,  averages  of  data  in  l*xl*  blocks  were  created  by  Levitus.  In 
order  to  overcome  biases  and  lack  of  data  and  to  improve  their 
representativeness,  a  smoothing  analysis  was  performed.  The  smoothing 
operation  was  a  weighted  average  operation  in  which  a  Gaussian  type  filter  was 
used  as  a  weighting  operator.  First,  approximate  values  for  each  grid  point 
were  assigned.  Then,  corrections  to  the  approximate  values  were  computed  as 
a  distance  weighted  mean  of  all  grid  points  within  the  area  specified  by  the 
radius  of  the  Gaussian  filter.  The  iterations  were  repeated  until  convergence. 
After  each  iteration,  the  resulting  data  set  was  smoothed  by  a  five  point 
smoother.  During  this  smoothing  operation  the  r'ldii  of  the  Gaussian  filter 
ranged  from  1540  km  to  770  km  resulting  in  the  total  elimination  of  any  signal 
with  wavelenghts  less  than  600  km,  while  wavelengths  between  600  km  and 
3000km  were  affected  with  changes  in  the  amplitudes  of  the  signal.  For 
example,  wavelengths  of  1000  km  had  a  reduction  in  amplitude  of  about  50% 
(Levitus,  1982,  Figure  11).  The  minimum  wavelengths  of  some  existing  signal  is 
of  the  order  of  800  km,  roughly  corresponding  to  an  harmonic  degree  20. 
Levitus  considers  that  the  resulting  large  scale  features  are  representative  of 
the  reed  ocean,  although  it  is  expected  that  some  local  differences  would  occur 
because  of  interannual  variability.  This  smoothed  data  set  was  used  to 
compute  the  annued  mean  anomaly  of  geopotential  thickness  of  the  0-2000  m 
layer  assuming  the  2000  m  level  as  a  level  of  no  motion  and  also  assuming 
geostrophy. 

The  Levitus  estimates  of  SST  were  available  to  us  in  l*xl*  means  in  the 
oceanic  areas  of  the  world.  The  global  mean  value  of  this  data  set  was  202  dyn 
cm.  First  analysis  of  that  data  indicated  the  existence  of  some  outliers  in  the 
west  Pacific  and  the  extremely  negative  values  of  SST  in  the  Mediterranean 
Sea.  Areas  with  no  data  were  found  in  the  North  Sea  and  along  many  of  the 
coastlines.  The  outliers  in  the  west  Pacific  were  substituted  by  the  global 
average.  For  the  Mediterranean  Sea,  estimates  from  a  map  of  Lisitzin  with 
respect  to  a  4000  meter  surface  (Lisitzin  1974,  p.l53)  were  used  to  substitute 
the  estimates  of  Levitus.  In  order  to  put  the  Lisitzin  estimates  in  the  same 
reference  system  as  the  Levitus  set,  their  mean  value  was  subtracted  and  the 
global  mean  value  of  Levitus  set  was  added. 

The  analysis  of  the  Levitus  SST  was  carried  out  similar  to  the  analysis  of 
the  Seasat-GEML2  derived  SST.  So,  a  new  mean  value  was  computed  (204  dyn 
cm)  and  subtracted  from  the  original  estimates  and  zeros  were  assigned  on 
land,  in  areas  with  no  data  and  polar  regions  (beyond  -70  and  70  degrees). 
Some  discontinuities  that  are  created  especially  in  areas  with  some  lack  of  data 
cannot  be  avoided,  although  they  are  not  important  since  they  mainly  distort 
the  higher  part  of  the  frequency  spectrum  and  so  they  do  not  significantly 
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affect  the  determination  of  the  large  scale  SST.  This  global  data  set  has  then 
been  harmonically  analysed,  using  procedures  appropriate  for  mean  values 
(Colombo,  1980).  The  harmonic  analysis  was  carried  out  up  to  degree  100, 
although  it  was  expected  that  beyond  degree  20  one  would  mainly  have  white 
noise.  The  power  spectrum  of  Levitus*  SST,  shown  in  Figure  Bl,  just  shows 
such  a  behavior.  Coefficients  up  to  degree  20  were  used  to  generate  the  SST 
map  shown  in  Figure  B4. 

Comparison  between  the  two  sets  up  to  degree  6  made  by  Engelis  (1985) 
indicates  that,  except  for  the  smoother  behavior  of  the  Levitus  set,  the  two 
sets  are  compatible  in  representing  qualitatively  the  major  features  of  the 
ocean  topography.  Indeed  only  their  amplitudes  differ,  while  there  is  no 
change  in  polarity,  with  the  exception  of  the  Indian  Ocean.  Moreover,  it  was 
seen  in  our  tests  that  contributions  of  degrees  7-20  in  the  Levitus  set  simply 
enhance  the  larger  scale  features  without  creating  any  changes  in  the  phase. 
These  facts  indicate  that  features  with  wavelengths  in  the  band  between  7  and 
20  harmonic  degrees  can  be  combined  with  the  lower  degree  coefficients 
derived  from  Seasat.  Tests  have  shown  that  such  an  addition,  although 
consistent,  is  not  very  significant  since  (mainly  because  of  the  smoothing  in 
the  Levitus  original  data)  the  total  power  in  the  band  7-20  is  334  cm*, 
accounting  only  for  10%  of  the  total  power  of  3340  cm*  of  the  Seasat  SST  up  to 
degree  6.  The  contribution  of  the  coefficients  between  7-20  degrees,  in  terms 
of  SST  itself,  has  a  maximum  effect  of  about  10cm  occuring  in  the  southern 
latitudes  while  in  most  of  the  oceanic  regions  it  hardly  exceeds  5  cm.  The 
power  spectrum  of  this  combined  set  is  given  in  Figure  B2  while  a  contour  map 
is  given  in  Figure  B5. 

For  the  gravity  anomaly  recovery  from  altimetry,  SST  can  be  used  in 
correcting  the  sea  surface  heights  before  their  usage  for  4g  recovery,  or 
alternatively  it  can  be  used  to  compute  the  effect  of  SST  coefficients  in  terms 
of  4g  correction  and  apply  that  correction  after  the  gravity  anomalies  are 
computed.  In  the  first  case  global  gridded  SST  estimates  have  to  be 
computed.  Then  by  one  simple  bilinear  interpolation  one  can  determine  the  SST 
correction  at  any  point  of  interest  and  apply  it  to  the  corresponding  sea 
surface  height  to  get  the  geoidal  undulation.  In  the  second  case  the  SST 
coefficients  are  used  to  compute  a  ig  correction,  using  the  traditional 
synthesis  for  gravity  anomalies.  Considering  that  the  second  method  is  much 
more  convenient  to  apply,  we  chose  it  for  the  correction  of  the  altimetrically 
derived  gravity  anomalies  due  to  SST  effects.  A  map  showing  the  magnitudes 
of  the  4g  correction  using  coefficients  up  to  n  =  20  is  shown  in  Figure  B6, 
while  the  power  spectrum  in  terms  of  4g  is  shown  in  Figure  B3. 
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Figure  B1 

Cumulative  and  By  Degree  Power  Spectrum  of  the  Levitus  SST.  Units  are  in  cm 
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Cumulative  and  By  Degree  Power  Spectrum  in  Terms  of  Anomaly  Correction  as  Computed  From  the 

Seasat— Levitus  SST  Coefficient  Set.  Units  are  in  mgal^. 


Surface  Topography  Based  on  the  Levitus  Data  Set  Using  a  Spherical  Harmonic  Expansion  to  Degree 

20  (Contour  Interval  is  20  cm) 


Figure  B5 

Topography  Based  on  the  Seasat/Levitus  Coefficients  to  Degree  '20  (Contour  Interval 


In  order  to  obtain  an  estimate  of  the  frequency  content  of  the  l/8*xl/8* 
gravity  amonalies  recovered  from  altimetry,  a  power  spectrum  analysis  of  eight 
profiles  is  performed.  These  profiles  are  60*  long  (481  data  points)  and  are 
taken  to  be  along  meridians,  so  that  there  is  a  uniform  spacing  in  all  of  them. 
The  locations  of  the  profiles  are  shown  in  Figure  Cl. 

Estimation  of  the  power  spectral  density  (or  simply  power  spectrum)  of 
discretely  sampled  deterministic  or  stochastic  signals  can  be  made  by  several 
methods  (Kay  et  al.,  1981).  One  approach  is  based  on  procedures  employing 
the  Fast  Fourier  Transform  (FFT)  technique  which  has  been  shown  to  give 
reasonable  results  for  a  large  class  of  signals  and  also  is  computationally 
efficient.  Problems  that  are  associated  with  this  approach  relate  to  the 
inability  to  distinguish  between  the  signal  and  the  noise  inherent  in  the 
sampled  data  and  to  the  low  resolution  of  the  recovered  spetrum. 

An  additional  problem  is  due  to  the  implicit  windowing  of  the  data  that 
occurs  when  processing  with  the  FFT.  Windowing  manifests  itself  as  leakage  in 
the  spectral  domain,  i.e.  energy  in  the  main  lobe  of  the  spectral  response  leaks 
into  the  sidelobes,  obscuring  and  distorting  other  spectral  responses  that  are 
present.  Appropriate  smoothing  of  the  spectrum  can  reduce  leakage  effects  at 
the  expense  of  resolution. 

Conventional  FFT  spectral  estimation  is  based  on  a  Fourier  series  model  of 
the  data.  Indeed  if  x(n)  is  a  series  of  N  data  at  N-1  equispaced  intervals  Ax, 
then  by  taking  its  Fourier  transform  we  obtain: 

N-1 

X(f)  =  x(n)e-J2"f‘nAx  (d) 

n=0 

where  f  is  the  frequency  (in  cycles  per  total  length  of  the  profile),  which 
according  to  the  well  known  sampling  theorem,  can  only  be  recovered  in  the 
interval  -(l/2Ax)<f^(l/2Ax).  (1/2ax  is  the  so-called  Nyquist  frequency).  Use  of 

the  FFT  permits  evaluation  of  the  above  equation  at  the  discrete  set  of  M 
equally  spaced  frequencies  fjii=roAf  for  m  =  0,  1,  ...M  and  Af  =  1/NAx.  Then  we 
obtain: 


Figure  Cl 

Location  of  Anomaly  Profiles  for  Which  Power  Spectra  Were  Computed 


Once  the  Fourier  transform  is  obtained,  the  power  spectrum  is  simply: 


I  Xml* 


m  =  0 


2 1  Xml*  0<io<M 


This  estimate  of  the  power  spectrum,  also  known  as  the  periodogram,  is 
usually  quite  erratic  due  to  leakage  and  has  to  be  smoothed.  Different 
smoothing  windows  can  be  applied.  One  of  them,  widely  suggested,  is  the 
Hanning  window  which  in  the  frequency  domain  is  simply  the  smoothing  of  the 
original  periodogram  as  follows: 


Pju  =  0.25P]||— I  +  O.SOPjn  +  0.25P]g-(-i 
Po  =  0.5Po  +  0.5Pi 
P^^  =  0.5Pm-i  0.5Pm 


An  alternative  procedure,  also  using  FFT  techniques  and  yielding  identical 
results  if  properly  applied,  is  the  one  via  an  autocovariance  function  estimate, 
also  known  as  the  Blackman  and  Tuckey  method  (Blackman  and  Tuckey,  1959). 
In  this  method  the  power  spectrum  is: 


ii 

Cxx(n)e-J2rrfn4x 


where  Cxx(n)  is  the  autocovariance  function  at  lags  nAx  up  to  a  maximum  of 
LAx  *  (N-l)Ax.  Again  f  can  only  be  recovered  in  the  interval 
-{1/2ax)  *  t  *  (1/2Ax).  Use  of  the  FFT  can  give  us  estimates  of  P(f)  at  the 
discrete  frequencies  fn,  =  mAf  where  Af  =  1/NAx,  up  to  f^  =  l/2Ax.  Then: 


This  method  yields  identical  results  with  the  previous  one  if  L  =  N-1.  The 
autocovariance  function  can  be  computed  by: 


N-n-1 


Cxx(n)  =  _  x(k+n)x(k) 


wmmw 


*-  %  .N  -%'l*. 


which  is  an  unbiased  estimate  of  the  true  autocovariance  function.  Since  the 
variance  of  this  estimate  tends  to  be  high  for  larger  lags  giving  an  erratic 
behaviour  to  Cxxi  is  usually  computed  up  to  1.  =  N/2  if  N  is  even  or  (N’-l)/2 
if  N  is  odd.  Additionally  a  smoothing  window  must  be  applied.  The  Hanning 
window  in  the  space  domain  is: 


W(n) 


ami 
+  cos-r I 


n  *  L 
otherwise 


(C9) 


To  smooth  the  autocovariance  function  the  window  weights  W{n)  are  multiplied 

The  Blackman  and  Tuckey  method  has  been  applied  to  the  eight  profiles 
under  investigation.  In  order  to  reduce  leakage  effects  due  to  existing  power 
at  wavelengths  larger  than  60*,  the  gravity  anomalies  defined  by  the  GEM9 
potential  coefficients  complete  up  to  degree  20  were  first  subtracted.  This 
operation  effectively  removed  all  constant  terms  and  substantially  reduced  any 
existing  trends. 

Since  the  spacing  (Ax)  between  the  gravity  anomalies  is  0M25  and  there 
are  481  data  samples,  the  frequency  increment  is: 


=  480^0. ^  cycle/60  degree 


(CIO) 


which  is  also  the  minimum  frequency  of  the  spectrum.  The  Nyquist  frequency 
is  4  cycles/degree.  On  the  other  hand,  since  the  wavelength  X  is  defined  as 
X  =  1/f,  the  maximum  wavelength  that  is  resolvable  is  X^g^j  =  60*  and  the 
minimum  is  Xg,j,i  =  0*.25. 

The  wavenumber  m  can  be  related  to  the  corresponding  spherical  harmonic 
degree  <,  using  the  following  ratio: 


<  mAf 

360* ■  1 


(Cll) 


or,  using  (CIO): 


I  =  6ffl 


(C12) 


In  principle  the  smallest  spherical  harmonic  degree  in  the  profiles  is  6  (i.e. 
m  =  1)  while  the  maximum  degree  is  1440  (i.e.  m  =  240).  Because  m  increases  in 
steps  of  1,  it  is  clear  from  (C12)  that  the  values  of  i  implied  by  the  profile  are 


determined  only  in  steps  of  6;  i.e.  the  frequency  resolution  in  terms  of 
harmonic  degrees  is  6.  If  a  reference  field  to  degree  20  is  usedi  the  minimum 
wavenumber  for  which  power  should  exist  in  the  spectrum  should  be  3.  If  a 
reference  field  to  degree  180  is  used,  the  minimum  wavenumber  is  30. 

After  the  power  spectra  for  all  the  profiles  have  been  computed  and 
plotted,  some  conclusions  can  be  drawn: 

First,  in  all  the  spectra  there  is  some  power  up  to  wavenumber  3.  This 
power  is  implied  by  the  residual  anomalies  with  respect  to  GEM9.  This 
contribution  is  expected  to  come  from  errors  in  the  GEM9  model,  errors  in  the 
altimetric  &g  due  to  long  wavelength  orbital  errors,  data  errors,  and  the 
neglect  of  sea  surface  topography. 

Second,  the  power  is  negligible  above  a  certain  wavenumber  in  the  range 
of  134-202  depending  on  the  profile.  This  indicates  that  the  0*125  data 
spacing  does  not  really  provide  us  information  up  to  the  corresponding 
maximum  wavenumber  which  is  240.  In  the  above,  negligible  contribution  is 
considered  (rather  arbitrarily)  to  be  less  than  0.3%  of  the  total  power.  The 
maximum  wavenumber  and  the  corresponding  wavelength  for  each  of  the 
profiles  are  shown  in  Table  Cl.  From  the  table  one  can  also  see  the  variances 
of  the  profiles  which  range  from  113  to  1521  mgal®. 

In  order  to  obtain  a  more  definitive  picture  of  the  power  spectrum  for  the 
eight  profiles,  their  average  spectrum  has  been  computed  and  shown  in  Figure 
C2.  Again  we  can  draw  similar  conclusions.  More  specifically,  the  power 
attributed  to  the  discrei>ancy  between  GEM9  and  altimetric  4g  still  exists  and 
amounts  to  about  77  mgal’  (cumulative  power  up  to  wavenumber  3).  The  formal 
accuracy  estimates  for  GEM9  give  a  total  error  variance  of  15  mgal^,  which  is 
not  consistent  with  the  discrepancy  of  77  mgal^  reported  above,  thus  indicating 
that  the  additioned  errors  (orbital  errors  etc)  are  expected  to  have  a 
significant  power,  of  the  order  of  60  rogaF. 


Table  Cl 

Power  Information  of  the  60*  Latitude  Profiles 


Profile 

Variance* 

ngaP 

Max  Frequency  of  non-negligible  power 

N.Lat 

Lon 

wavenumber 

wavelength 

61 

175 

853 

160 

0:375 

61 

185 

1013 

134 

0:447 

61 

335 

393 

137 

0:437 

-1 

211 

201 

0:300 

-1 

90 

208 

176 

0:340 

-1 

306 

171 

0:350 

-1 

113 

202 

0:300 

-1 

335 

1521 

136 

0:440 

Average 

577 

155 

0:380 

*  with  respect  to  the  GEM9  gravity  field  to  degree  20 
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Figure  C2 

Average  Power  Spectrum  From  Eight  60*  Latitude  Profiles. 

(Units  of  frequency  are  cycles/(J0  degrees) 


In  a  similar  way,  it  is  observed  that  about  75X  of  the  total  power  is 
contained  in  frequencies  corresponding  to  degrees  up  to  <inax  - 
cycles/degree)  and  about  90%  to  frequencies  corresponding  to  tmax  =  360  (1 
cycle/degree).  In  the  average  power  spectrum  the  contribution  of  power  is 
negligible  (less  than  0.3%)  above  the  wavenumber  155  (2.6  cycles/degree)  or 
above  harmonic  degree  930.  This  corresponds  to  a  minimum  wavelength  of 
about  0*.38  or  in  a  data  spacing  of  0*.19.  This  minimum  wavelength  is 
consistent  with  the  one  found  by  Brammer  and  Sailor  (1982).  From  a  coherent 
analysis  of  overlapping  passes  of  SEASAT  they  have  established  that  the 
minimum  wavelength  that  can  be  recovered  from  altimetry  is  of  the  order  of  40 
km  (0*.36).  Smaller  wavelengths  cannot  be  recovered  since  the  signal  to  noise 
ratio  is  unfavorable  and  so  the  signal  is  completely  contaminated  by  noise.  In 
our  spectrum  estimation  the  power  at  these  frequencies  is  zero  because  the 
noise  and  the  contaminated  signal  have  been  effectively  removed  by  the  least 
squares  collocation. 

Considering  the  above  discussion  the  selected  spacing  of  0M25  is  justified 
since  a  wider  spacing,  say  0*.25,  would  result  in  some  loss  of  information.  This 
argument  becomes  stronger  if  we  consider  that  the  aliasing  effects,  which 
always  arise  during  the  discretization  of  a  continuous  non-bandlimited  signal, 
would  be  larger  if  a  wider  spacing  had  been  chosen. 

The  above  discussions  yield  the  conclusion  that  these  power  spectral 
estimates  are  qualitatively  correct.  It  is  believed  that  analysis  of  more  profiles 
will  give  more  reliable  results  although  they  are  not  expected  to  considerably 
differ  from  the  ones  analysed  above. 
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